Index of /rozprawy2/10608

Pełen tekst

(1)Akademia Górniczo-Hutnicza im. Stanisława Staszica. Wydział Energetyki i Paliw Katedra Podstawowych Problemów Energetyki. Rozprawa doktorska. ANALIZA WPŁYWU POLA MAGNETYCZNEGO NA OBSZAR. PRZEJŚCIOWY. POMIĘDZY. LAMINARNĄ A TURBULENTNĄ. mgr inż. Łukasz Pyrda. Promotor: Prof. dr hab. inż. Janusz S. Szmyd. Kraków 2013. KONWEKCJĄ.

(2) II.

(3) Mojej rodzinie. III.

(4) IV.

(5) Spis treści Spis symboli .......................................................................................................... VII 1 Wstęp .............................................................................................................................1 1.1 Właściwości magnetyczne materiałów ...................................................................1 1.2 Konwekcja naturalna ...............................................................................................3 1.3 Konwekcja termomagnetyczna płynów ..................................................................7 1.4 Dotychczasowy stan badań ...................................................................................10 1.5 Cel pracy ...............................................................................................................13 2 Model matematyczny konwekcji termomagnetycznej ................................................15 2.1 Magnetyzacja materiałów paramagnetycznych ....................................................15 2.2 Siły masowe ..........................................................................................................17 2.3 Równania warunków .............................................................................................19 3 Analiza numeryczna ....................................................................................................21 3.1 Modelowanie numeryczne ....................................................................................21 3.2 Schematy numeryczne ..........................................................................................22 3.3 Właściwości płynu ................................................................................................23 3.4 Rozkład pola magnetycznego ...............................................................................24 3.5 Geometria i warunki brzegowe .............................................................................27 3.6 Wprowadzenie do wyników obliczeń ...................................................................28 3.7 Wariant obliczeń nr 1 (Pr=20)...............................................................................29 3.8 Wariant obliczeń nr 2 (Pr=100).............................................................................41 3.9 Wariant obliczeń nr 3 (Pr=146).............................................................................53 3.10 Wariant obliczeń nr 4 (Pr=584)...........................................................................65 3.11 Analiza porównawcza wyników numerycznych .................................................77 4 Analiza eksperymentalna .............................................................................................81 4.1 Wprowadzenie.......................................................................................................81 4.2 Układ eksperymentalny.........................................................................................81 V.

(6) 4.3 Płyny robocze ........................................................................................................83 4.4 Procedura pomiarowa............................................................................................88 4.5 Pomiary termiczne ................................................................................................89 4.6 Wprowadzenie do wyników eksperymentalnych..................................................91 4.7 Analiza eksperymentalna płynu P0 (∆θ=5 [K]) ....................................................93 4.8 Analiza eksperymentalna płynu P50 (∆θ=5 [K]) ................................................102 4.9 Analiza eksperymentalna płynu P50 (∆θ=11 [K]) ..............................................111 4.10 Analiza eksperymentalna płynu P80 (∆θ=5 [K]) ..............................................120 4.11 Analiza porównawcza wyników eksperymentalnych .......................................128 5 Wnioski ogólne ..........................................................................................................135 Bibliografia ........................................................................................................... 137 Dodatek A ............................................................................................................. 143 Dodatek B ............................................................................................................. 147 Dodatek C ............................................................................................................. 157. VI.

(7) Spis symboli A. powierzchnia [m2]. aloss. współczynnik proporcjonalności w równaniu określającym straty ciepła [W/K]. b. indukcja magnetyczna [T]. b0. indukcja magnetyczna pola zewnętrznego [T]. b0 max. indukcja magnetyczna w centrum magnesu nadprzewodzącego [T]. ℂ. liczby zespolone [-]. C. stała Curie [K·m3/kg]. Cbal. stała kalibracyjna wagi magnetycznej [-]. Cp. ciepło właściwe [J/kg·K]. D. wymiar charakterystyczny układu [m]. ds. wektor elementu przewodnika prądem [m]. e. liczba Eulera, e = 2,718... [-]. fB. grawitacyjna siła wyporu [N/m3]. fM. magnetyczna siła wyporu [N/m3]. f gr. siły masowe wywołane działaniem pola grawitacyjnego [N/m3]. f mg. siły masowe wywołane działaniem pola magnetycznego [N/m3]. f. częstotliwość [Hz]. g. przyśpieszenie ziemskie [m/s2]. H. natężenie pola magnetycznego [A/m]. h/d. współczynnik kształtu, wysokość do szerokości [-]. h1. grubość ścianki horyzontalnej naczynia pomiarowego [m]. h2. grubość ścianki adiabatycznej naczynia pomiarowego [m] VII.

(8) i. jednostka urojona [-]. I. natężenie prądu [A]. k. numer harmonicznej [-]. Ls. wysokość próbki w zestandaryzowanej próbówce [m]. M. magnetyzacja, namagnesowanie [A/m]. M max. namagnesowanie nasycenia [A/m]. mS. masa próbki w próbówce [kg]. N. liczba próbek sygnału [-]. n. numer próbki sygnału [-]. p. ciśnienie [N/m2]. ℝ. liczby rzeczywiste [-]. RS. odczyt wagi magnetycznej dla próbówki z próbką [-]. RS 0. odczyt wagi magnetycznej dla pustej próbówki [-]. RT. różnica odczytów wagi magnetycznej dla próbki i pustej próbówki [-]. r. wektor wodzący [m]. r. długość wektora wodzącego [m]. Q. strumień ciepła [W]. Qcond. strumień ciepła wymienianego na drodze przewodzenia [W]. Qconv. strumień ciepła wymienianego na drodze konwekcji [W]. Qheater. strumień ciepła dostarczony do układu [W]. Qloss. całkowity strumień ciepła oddawany do otoczenia w postaci strat [W]. Qnet _ cond. strumień ciepła wymienianego na drodze przewodzenia netto [W]. Q net _ conv. strumień ciepła wymienianego na drodze konwekcji netto [W]. Qtheor _ cond. teoretyczny strumień ciepła wymienianego na drodze przewodzenia [W]. VIII.

(9) t. czas [s]. U. napięcie elektryczne [V]. u. prędkość [m/s]. xc. długość naczynia pomiarowego [m]. yc. szerokość naczynia pomiarowego [m]. z. odległość od centrum magnesu nadprzewodzącego [m]. zc. wysokość naczynia pomiarowego [m]. Symbole greckie α. współczynnik dyfuzyjności termicznej [m2/s]. β. współczynnik rozszerzalności termicznej [1/K]. ε. współczynnik dyssypacji lepkiej [m2/s3]. θ. temperatura [K]. θ0. temperatura referencyjna [K]. θc. temperatura ścianki chłodnej [K]. θ Curie. temperatura Curie [K]. θh. temperatura ścianki ciepłej [K]. θ∞. średnia temperatura otoczenia układu [K]. λ. współczynnik przewodzenia ciepła [W/m·K]. η. skala długości Kolmogorova [m]. η A ,η B ,η C. współczynniki. wielomianu. opisującego. zależność. dynamicznego. współczynnika lepkości od temperatury [-]. µ. dynamiczny współczynnik lepkości [kg/m·s]. µθ. dynamiczny współczynnik lepkości (zależny od temperatury) [kg/m·s]. µ0. przenikalność magnetyczna próżni, µ 0 = 4π ⋅ 10 −7 [H/m] IX.

(10) µm. przenikalność magnetyczna [H/m]. τη. skala czasu Kolmogorova [s]. ν. kinematyczny współczynnik lepkości [m2/s]. π. liczba Pi, π = 3,14... , [-]. ρ. gęstość [kg/m3]. ρ0. gęstość w temperaturze θ 0 [kg/m3]. χ. podatność magnetyczna [-]. χm. masowa podatność magnetyczna [m3/kg]. χm0. masowa podatność magnetyczna w temperaturze θ 0 [m3/kg]. Liczby podobieństwa γ. χ b0 max parametr gamma, liczba opisująca magnetyzację płynu [-], γ = µm gD. Gr. liczba Grashofa [-], Gr=. Nu. liczba Nusselta [-], Nu=. Qconv Qcond. Pr. liczba Prandtla [-], Pr =. ν α. Ra T. termiczna liczba Rayleigha [-], Ra T =. Ra M.  1   γ gβ (θ h − θc ) D3  magnetyczna liczba Rayleigha [-], Ra M = 1 +   2αν   βθ0  . 2. gβ (θ h − θ c ) D3. ν2. gβ (θ h − θ c ) D 3. αν. Ra TM termomagnetyczna liczba Rayleigha [-], Ra TM = =Ra T +Ra M. X.

(11) 1 Wstęp. 1 Wstęp 1.1 Właściwości magnetyczne materiałów Equation Chapter 1 Section 1 Umieszczenie materiału w polu magnetycznym powoduje, że siła magnetyczna zaczyna działać na znajdujące się w nim elektrony. Materiały mogą reagować na pole magnetyczne w różnorodny sposób, w zależności od wielu czynników, np. od struktury atomowej czy molekularnej materiału. Stopień namagnesowania materiału w zewnętrznym polu magnetycznym jest nazywany podatnością magnetyczną ( χ ), zdefiniowaną równaniem [1]:. χ=. M , H. (1.1). gdzie M jest namagnesowaniem materiału, a H jest natężeniem pola magnetycznego. Masowa podatność magnetyczna χm jest ilorazem podatności magnetycznej i gęstości ( ρ ) [1] :. χm =. χ . ρ. (1.2). Indukcja magnetyczna b zależy od natężeniem pola magnetycznego H oraz od przenikalności magnetycznej danego materiału, co można zdefiniować następująco: b = µmH. (1.3). przy relacji:. µm = µ0 (1 + χ ) .. (1.4). gdzie: µm to przenikalność magnetyczna, natomiast µ0 to przenikalność magnetyczna próżni. W zależności od reakcji materiałów na przyłożone zewnętrzne pole magnetyczne, można je podzielić na trzy najważniejsze grupy: diamagnetyki,. 1.

(12) paramagnetyki i ferromagnetyki. Taki podział pierwiastków został przedstawiony na rysunku 1−1 [2].. Rys. 1-1 Podział pierwiastków ze względu na właściwości magnetyczne w temperaturze 300 [K] [2]. Diamagnetyzm jest cechą wszystkich materiałów, jednak jest to efekt tak słaby, że jest „maskowany” przez wypadkowy moment magnetyczny atomów występujący w paramagnetykach umieszczonych. i. ferromagnetykach. w zewnętrznym. polu. [3].. W. materiałach. magnetycznym,. diamagnetycznych. generowane jest. pole. magnetyczne skierowane przeciwnie do tego pola, niezależnie od temperatury. Dla tych materiałów podatność magnetyczna ośrodka jest ujemna i mniejsza od jedności [4]. Diamagnetyki samorzutnie nie wykazują właściwości magnetycznych. Są odpychane przez magnes, tak, że w przypadku przyłożenia bardzo silnego zewnętrznego pola magnetycznego mogą lewitować. Paramagnetyzm przejawiają materiały zawierające pierwiastki ziem rzadkich, pierwiastki przejściowe i aktynowce. Po umieszczeniu w zewnętrznym polu magnetycznym, elementarne dipole atomowe takich materiałów będą usiłowały ustawić się w kierunku zgodnym z kierunkiem pola, co nazywane jest paramagnetyzmem [3]. Tworzą one wewnętrzne pole magnetyczne o kierunku zgodnym z kierunkiem. 2.

(13) 1 Wstęp przyłożonego pola magnetycznego. Materiały te posiadają dodatnią podatność magnetyczną, mniejszą od jedności [4]. Są lekko przyciągane przez źródło pola magnetycznego. Ferromagnetyzm jest własnością żelaza, niklu i kobaltu oraz stopów tych pierwiastków [3]. Materiały ferromagnetyczne posiadają bardzo dużą dodatnią podatność magnetyczną. Wykazują bardzo silne przyciąganie do zewnętrznego pola magnetycznego i potrafią zachować własności magnetyczne po jego usunięciu zewnętrznego pola magnetycznego, wtedy otrzymuje się tzw. magnesy trwałe [4]. Ich silne właściwości magnetyczne istnieją dzięki obecności domen magnetycznych, w których duże ilości atomów posiadają spinowe momenty magnetyczne elektronów ułożone równolegle. W przypadku, gdy temperatura ferromagnetyku wzrośnie powyżej krytycznej wartości, zwanej temperaturą Curie to materiał staje się paramagnetykiem.. 1.2 Konwekcja naturalna. Konwekcja naturalna stanowi taki sposób wymiany ciepła, w którym ruch płynu wywołany jest siłami masowymi, najczęściej siłami ciężkości, działającymi na części płynu o różnych gęstościach. W otoczeniu ciała wymieniającego ciepło występują różnice temperatury, powodujące różnicę gęstości otaczającego płynu, która z kolei jest przyczyną powstawania siły wyporu. Konwekcja naturalna płynu w zamkniętej przestrzeni jest tak zróżnicowana, jak różne są rozpatrywane geometrie i ich orientacja w przestrzeni. Rozpatrując to zjawisko w zależności od możliwego zastosowania, można je podzielić na dwie duże klasy: geometrie z ogrzewaną dolną ścianką i chłodzoną górną oraz pozostałymi adiabatycznymi, geometrie z ogrzewaną ścianką boczną i drugą przeciwległą chłodzoną oraz pozostałymi adiabatycznymi. W celu identyfikacji możliwego typu konwekcji stosuje się bezwymiarowe liczby podobieństwa: liczbę Grashofa, opisującą stosunek siły wyporu do sił lepkości danego płynu:. Gr =. gβ (θ h − θc ) D 3. ν2. ,. (1.5). 3.

(14) gdzie:. D. - wymiar charakterystyczny układu, g. - przyśpieszenie ziemskie,. ν - kinematyczny współczynnik lepkości, θ h - temperatura ścianki ciepłej, θc - temperatura ścianki zimnej, β - współczynnik rozszerzalności termicznej, liczbę Prandtla, opisującą stosunek dyfuzyjności pędu do dyfuzyjności termicznej:. Pr = gdzie:. Cp. -. ciepło. właściwe,. µ. Cp µ. λ -. =. ν , α. (1.6). dynamiczny. współczynnik. lepkości,. λ - współczynnik przewodzenia ciepła, α - współczynnik dyfuzyjności termicznej, oraz (termiczną) liczbę Rayleigha będącego iloczynem liczby Grashofa i Prandtla:. Ra T = Gr Pr =. gβ (θ h − θ c ) D 3. αν. .. (1.7). W celu opisu intensywności wymiany ciepła w danym układzie, stosuję się bezwymiarową liczbę Nusselta, która wyraża stosunek wymiany ciepła na drodze konwekcji ( Qconv )do teoretycznej wymiany ciepła na drodze przewodzenia ( Qcond ): Nu=. Qconv . Qcond. (1.8). Na podstawie znanej liczby Prandtla oraz Rayleigha, można aproksymować liczbę Nusselta korzystając z korelacji otrzymanych eksperymentalnie przez wielu badaczy. Przykładowe korelacje liczby Nusselta zestawiono w tabeli 1-1. Należy podkreślić, iż przedstawione korelacje w tabeli 1-1 nie są uniwersalne, gdyż mają zastosowanie jedynie dla wybranego rodzaju geometrii oraz właściwości fizycznych płynu. Rozważając geometrię sześcienną oraz umieszczony w niej płyn o liniowej zależności gęstości od temperatury, można zaobserwować, w zależności od rozpatrywanych warunków brzegowych różne efekty. W konfiguracji z ogrzewaną. ścianką dolną i chłodzoną górną (tzw. konfiguracja Rayleigha- Bénarda) rozpatrywany układ jest niestabilny i konwekcja laminarna pojawia się po przekroczeniu krytycznej. 4.

(15) 1 Wstęp liczby Rayleigha ( Ra T >1708 [5]). Sytuację taką przedstawiono na rysunku 1−2a i 1−2b, gdzie niebieskim kolorem oznaczono płyn chłodny, natomiast czerwonym ciepły. Tabela 1-1 Wybrane, eksperymentalne korelacje liczby Nusselta względem liczby Rayleigha, Grashofa oraz Prandtla. Autorzy. Korelacja. Przedział Ra. Pr. Chu i Goldstein [6]. Nu = 0,183Ra 0,278. 2 ⋅ 10 5 − 10 8. 6. Nu = 0,14Ra 0,26. 4 ⋅10 6 − 4 ⋅ 108. 0, 025. Nu = 0, 44Ra 0,20. 4 ⋅108 − 2 ⋅10 9. 0, 025. Nu = 0,145Ra 0,292. 3 ⋅10 − 7 ⋅10. 6,5. Nu = 0, 069Ra1/3 Pr 0,074. 3 ⋅10 5 − 7 ⋅10 9. 0, 02 − 11560. Nu = 1 + 0, 096Ra1/3. 3 ⋅10 5 − 4 ⋅10 7. 0, 6 − 1, 4. Nu = 1 + 0, 2Ra 282. 4 ⋅ 10 7 − 1011. 0, 6 − 1, 4. Nu = 0,124Ra 0,309. 10 6 − 1017. 0, 7 − 30. Nu = 0,147Ra 0,247. 2 ⋅ 10 4 − 5 ⋅10 5. 0, 024. Nu = 0,131Ra 0,3. 3 ⋅10 4 − 3 ⋅ 10 6. 6,8. Nu = 0,184Ra 0,281. 4 ⋅ 10 3 − 3 ⋅10 6. 200. Silveston [12]. Nu = 0,1Gr 0,3 ⋅ Pr 0,36. 4, 4 ⋅10 4 − 107. 0, 02 − 3000. Nu = 0,196Ra. 5,5. Sommerscales [13]. 3 ⋅10 7 − 3 ⋅108 4 ⋅ 10 5 − 10 8. 18. Cioni [7]. Dropkin [8] Heslot [9] Niemela[10]. Rosby [11]. 0,283. 8. 9. W przypadku konfiguracji z ogrzewaną ścianką boczną i przeciwległą chłodzoną, układ jest zawsze niestabilny i konwekcja laminarna występuje w przypadku różnej temperatury pomiędzy ściankami bocznymi (rys. 1−2c). Ostatnim przedstawionym układem na rysunku 1−2d jest układ zawsze stabilny, w którym niezależnie od różnicy temperatury występuje jedynie przewodzenie.. 5.

(16) Rys. 1-2 Przykładowe mapy temperatury w geometrii sześciennej: a) i b) układ niestabilny po przekoczeniu czeniu krytycznej liczby Rayleigha, c) układ zawsze niestabilny, d) układ stabilny. Wraz ze wzrostem liczby Rayleigha zmienia się równieżż charakter wymiany ciepła, która z konwekcji laminarnej przechodzi w turbulentną, turbulentną, przechodząc przechodz poprzez obszar przejściowyy noszący noszący cechy zarówno konwekcji laminarnej jak i turbulentnej. W obszarze konwekcji turbulentnej następuje nast puje intensyfikacja wymiany ciepła oraz mieszania płynu. Typ konwekcji zależy zale od liczby Prandtla a także że geometrii wybranego układu. Na rysunku 1− −3 został przedstawiony schemat obrazujący ący obszary przejścia przej pomiędzy konwekcjąą naturalną naturaln a turbulentną dla (h/d ≪ 1) klasycznego eksperymentu Rayleigha-Bénarda. Gdy liczba Rayleigha jest o rzędy wielkości większa ksza od krytycznej wartości, warto konwekcja onwekcja przejawia charakter turbulentny. Rdzeń Rdze płynu posiada praktycznie temperaturę średnią (θ h + θ c ) / 2 , podczas gdy spadek temperatury występuje wyst w cienkich, horyzontalnych warstwach przyściennych przy ciennych i jest równy. (θ h − θ c ) / 2 .. Turbulencja urbulencja w warstwach horyzontalnych jest wywołana wtrąceniami wtrąceniami meandrujących, meandruj termicznych strumieni płynu, zarówno wyrastających wyrastaj cych z podgrzewanej warstwy dolnej, jak i spadających cych z chłodzonej warstwy górnej [14].. 6.

(17) 1 Wstęp. Rys. 1-3 Trójwymiarowa reprezentacja zależności Ra, Pr i proporcji geometrycznych (wysokość/szerokość) boków. Rezultat na płaszczyźnie Ra - Pr jest wyznaczony dla (h/d ≪ 1) klasycznego przypadku Rayleigha-Bénarda [15]. Turbulencja pojawia się, gdy w przepływającym płynie, siły wyporu, bezwładności lub inne przezwyciężą lepkość płynu, która przeciwdziała zaburzeniom i niestabilnościom przepływu. Jeśli bezwładność płynu i siła masowa są wystarczająco duże w porównaniu do lepkości płynu, to wzmacniają losowe zaburzenia, występujące w każdym przepływie. Zaburzenia te rosną, stają się niestabilne i nieliniowe, oddziałują na siebie wzajemnie i skupiają się w chaotyczny ruch.. 1.3 Konwekcja termomagnetyczna płynów. Konwekcja termomagnetyczna płynu nieprzewodzącego prądu elektrycznego, umieszczonego w niejednorodnym polu magnetycznym, na który działają siły ciężkości, wynika z sumy sił oddziaływania grawitacyjnego i magnetycznego. W zależności od rozpatrywanego materiału magnetycznego i ułożenia względem siebie wektorów siły magnetycznej i grawitacyjnej można otrzymać różne efekty wyporu grawitacyjno-magnetycznego. 7.

(18) Wpływ niejednorodnego pola magnetycznego na konwekcję termiczną płynów zależy od zmiany ich własności magnetycznych w funkcji temperatury. W zależności od reakcji materiałów na przyłożone zewnętrzne pole magnetyczne, można je podzielić na trzy główne grupy [16]: a) diamagnetyki – podatność masowa nie zmienia się z temperaturą, pole magnetyczne ma znikomy wpływ na konwekcję termiczną cieczy, b) paramagnetyki – podatność magnetyczna jest odwrotnie proporcjonalna do temperatury, pole magnetyczne wywiera silny wpływ na konwekcję termiczną cieczy i gazów, c) ferromagnetyki – podatność magnetyczna jest funkcją nieliniową temperatury poniżej temperatury Curie oraz nieróżniczkowalną względem temperatury, powyżej temperatury Curie, gdzie materiał staje się paramagnetykiem i jego podatność magnetyczną opisuje prawo Curie-Weissa [17]. Ferromagnetyki są najczęściej płynami przewodzącymi prąd elektryczny. W takich płynach pole magnetyczne indukuje potencjały elektryczne. Rozpatrywana jest dodatkowo siła Lorenza, która jest znacznie silniejsza od siły wyporu magnetycznego (zwanego również siłą Kelvina). Badaniem wpływu pola magnetycznego na płyny przewodzące prąd elektryczny zajmuje się magnetohydrodynamika (MHD) [18]. Schematyczną zależność masowej podatności magnetycznej od temperatury dla różnych typów substancji magnetycznych przedstawiono na rysunku 1−4. Na rysunku 1−5 przedstawiono dwie przykładowe konfiguracje sił wyporu: magnetycznego i grawitacyjnego. Kolor niebieski oznacza chłodny płyn, czerwony ciepły w odniesieniu do temperatury referencyjnej, natomiast elipsą, oznaczona została pozycja centrum pola magnetycznego. W przypadku przedstawionym na rysunku 1−5a kierunki i zwroty działania sił wyporu są zgodne, czego efektem jest intensyfikacja wymiany masy i ciepła. W przykładzie zaprezentowanym na rysunku 1−5b wektory tych sił są skierowane przeciwnie, co przy odpowiednio dużej sile wyporu magnetycznego umożliwia zrównoważenie grawitacyjnej siły wyporu lub nawet jej przezwyciężenie.. 8.

(19) 1 Wstęp. Rys. 1-4 Schematyczna zależność masowej podatności magnetycznej od temperatury dla: a) diamagnetyków, b) paramagnetyków (zgodnie z prawem Curie), c) ferromagnetyków [16]. Rys. 1-5 Konfiguracja siły wyporu: magnetycznego i grawitacyjnego: a) zgodny kierunku i zwrot sił, b) zgodny kierunek i przeciwny zwrot sił. 9.

(20) 1.4 Dotychczasowy stan badań. Konwekcja naturalna Pierwsze ilościowe eksperymenty związane z laminarną konwekcją naturalną przeprowadził Henri Bénard [19] w 1900 roku. Dotyczyły one konwekcji w przestrzeni pomiędzy dwoma równoległymi płytami o nieskończonych wymiarach w stosunku do pozostałych powierzchni. Dolna płyta była ogrzewana. Bénard w przeprowadzonym eksperymencie zaobserwował pojawienie się komórek heksagonalnych. W 1916 roku Lord Rayleigh [20] opublikował pracę odnoszącą się do eksperymentu Bénarda, w której sformułował podstawy teoretyczne zjawiska znanego dzisiaj jako konwekcja Rayleigha-Bénarda. Zdefiniował również parametr bezwymiarowy (później nazwany liczbą Rayleigha), który aby zaszła konwekcja, musi przekroczyć pewną krytyczną wartość. W kolejnych latach prace eksperymentalne i teoretyczne nad tym zagadnieniem w konfiguracjach sześciennych, prowadzone były przez wielu naukowców m. in. Jeffreys [21], Schmidt [22], Pellew [23], Silveston [12], Dropkin [8], Sommerscales [13], Rosby [11], Busse [24], Krishnamurti [25]. Badania wymienionych autorów nad zjawiskiem konwekcji koncentrowały się m. in na identyfikacji typu konwekcji, jej wizualizacji jakościowej i ilościowej oraz próbie wyznaczenia korelacji liczb bezwymiarowych na podstawie otrzymanych danych eksperymentalnych. Większość badań przeprowadzonych do tego okresu opisuje zachowanie płynu, jako turbulentne, gdy tylko przepływ staje się nieokresowy. Aktualnie wiadomo, że jest to stan znacząco różny od rozwiniętej konwekcji turbulentnej. W 1987 Heslot [9] zaproponował klasyfikację tego obszaru na trzy regiony: obszar przejściowy, lekką oraz silną turbulencję. Region lekkiej turbulencji w przedziale 10 5 < Ra < 4 ⋅ 10 7 , natomiast silnej w przedziale 4 ⋅10 7 < Ra < 6 ⋅1011 . Podział ten zaproponował na podstawie analizy statystycznej otrzymanych wyników eksperymentalnych. Castaing [26] przeprowadził badania eksperymentalne przy wykorzystaniu helu, utrzymywanego w temperaturze 5 [K] oraz układu badawczego o stosunku ścianek równemu jedności (h/d=1). Badania wykonano w rejonie silnej turbulencji ( 4 ⋅ 10 7 <Ra < 1012 oraz 0, 65 < Pr < 1, 5 ). Zaobserwowano jednorodne zachowanie płynu, w całym badanym. obszarze, jednakże w przeciwieństwie do wcześniejszych badań, korelacja liczby Rayleigha z liczbą Nusselta była bliższa Ra 2/ 7 niż Ra 1/3 . Najnowsze badania 10.

(21) 1 Wstęp przeprowadzane przez Niemela [10] i innych, koncentrują się na obszarze silnej turbulencji ( 10 6 < Ra < 1017 ). Badania te również przeprowadzono przy wykorzystaniu helu, utrzymywanego w temperaturze 5 [K]. Rezultatem tych badań jest m. in. korelacja Nu = 0,124Ra 0,309 , która nie potwierdza wcześniejszych teorii oraz korelacji.. Podstawowe definicje i charakterystyki związane z konwekcją naturalną zawarte są w opracowaniach książkowych np. Jaluria [27], Bejan [14], Chandrasekhar [28], Drazin i Reid [29]. Pierwsze symulacje numeryczne laminarnej konwekcji naturalnej w konfiguracji Rayleigha-Bénarda przeprowadził Deardoff [30] w 1964 roku. Przed uruchomieniem symulacji wprowadzano niewielkie zaburzenia do liniowego gradientu temperatury po czym rozpoczynano obliczenia. Po pewnym czasie pole temperatury ulegało stabilizacji, a otrzymane wyniki porównywano z wynikami eksperymentalnymi. Grötzbach [31] w 1984 roku przedstawił wyniki bezpośredniej symulacji numerycznej (DNS), zarówno dla obszaru laminarnego jak i turbulentnego ( Ra=3,8 ⋅105 ). Symulację wykonał na stosunkowo rzadkiej siatce obliczeniowej 32x32x16, a rezultaty porównał z danymi eksperymentalnymi. Stwierdził również, że należy zagęścić siatkę obliczeniową w rejonie przyściennym w celu dokładniejszego zbadania rozkładu temperatury. Henkes [32] wykazał, że obszar konwekcji przejściowej jest trudny do symulacji numerycznych, z powodu braku odpowiednich modeli turbulencji, co wymaga zastosowania dokładnych siatek numerycznych oraz odpowiednio małego kroku obliczeniowego. Janssen [33] przeprowadził szereg obliczeń numerycznych konwekcji laminarnej, przejściowej i turbulentnej i przedstawił wpływ liczby Prandtla na mechanizm niestabilności oraz przejście pomiędzy typami konwekcji. Kerr [34] przedstawił serię wyników obliczeń DNS wykonanych m.in. dla powietrza przy Ra = 2 ⋅ 10 7 , gdzie zaobserwowano silną turbulencję w przedziale 5 ⋅ 10 4 ≤ Ra ≤ 2 ⋅10 7 .. W publikacji [34] poruszono również szereg problemów związanych z modelowaniem konfiguracji Rayleigha-Bénarda, m. in. wpływ różnych warunków brzegowych. Ze względu na to, że obliczenia DNS nie są wystarczająco wydajne, aby stosować sztywne ścianki na więcej niż jednej osi zastosowano geometrię o stosunku boków 6:6:1. Do obliczeń wykorzystano siatkę złożoną z 288x288x96 węzłów. Kenjeres [35] określił limit symulacji DNS do Ra ≈ 1010 , ze względu na wymaganą ilość węzłów obliczeniowych niezbędną, aby symulować wszystkie skale turbulencji. Dlatego,. 11.

(22) zamiast niej zastosował symulację TRANS (transient RANS), która wykorzystuje uśrednione równania Naviera-Stokesa. Dzięki temu przeprowadził symulację dla 10 6 ≤ Ra ≤ 2 ⋅ 1014 , której wyniki porównano z dostępnymi wynikami symulacji DNS. oraz danymi eksperymentalnymi.. Magnetyzm i konwekcja magnetyczna Pierwsze badania potwierdzające związek pomiędzy zjawiskami elektrycznymi i magnetycznymi przeprowadził w 1820 roku Oersted. Odkrył, że wokół prostoliniowego przewodnika, przez który przepływa prąd elektryczny, występuje pole magnetyczne. W tym samym roku Jean-Baptise Biot i Felix Savart [36] wykazali, że indukcja magnetyczna w danym punkcie jest wprost proporcjonalna do natężenia prądu elektrycznego płynącego w przewodniku oraz zależy od kształtu i rozmiarów przewodnika, w którym płynie prąd. W 1847 roku Faraday [37] wykazał, że wszystkie substancje posiadają własności magnetyczne. Przeprowadził szereg eksperymentów z różnymi substancjami, przechodząc od badań wpływu pola magnetycznego na płomień do obserwacji wpływu pola magnetycznego na gazy takie jak tlen czy wodór. Trzy lata później, Lord Kelvin opisał przenikalność magnetyczną i podatność magnetyczną i ich związek z natężeniem pola magnetycznego. W 1895 roku Pierre Curie [38] w swojej pracy doktorskiej wykazał, że podatność magnetyczna paramagnetyków jest odwrotnie proporcjonalna do temperatury oraz przedstawił brak wpływu temperatury na podatność magnetyczną diamagnetyków. Wynalezienie nadprzewodników działających powyżej temperatury wrzenia ciekłego azotu (77 [K]) w 1987 roku, umożliwiło konstrukcję magnesów nadprzewodnikowych [39]. W 1991 roku Braithwaite wraz z zespołem [40] przedstawił wpływ silnego niejednorodnego pola magnetycznego na proces konwekcji naturalnej płynów paramagnetycznych oraz opis matematyczny sił działających na taki płyn, w którym występują różnice temperatur. W 1998 roku, Huang i Edwards [41] opisał możliwość intensyfikacji i spowolnienia konwekcyjnego procesu wymiany ciepła w nieprzewodzącym prądu elektrycznego płynie paramagnetycznym. Duży wkład w badania nad konwekcją termomagnetyczną włożył Ozoe wraz z zespołem, a wyniki wieloletnich prac opublikował w książce „Magnetic Convection” [42]. Tagawa i in. [43] przedstawili model matematyczny oraz wyniki obliczeń numerycznych dla konwekcji termomagnetycznej gazu paramagnetycznego w układzie sześciennym. 12.

(23) 1 Wstęp Tagawa i in. [44] wykonali również podobne obliczenia wykonane dla cieczy diamagnetycznej. Bednarz i in. [45], [46], [47], [48], [49] wykonali szereg badań nad wpływem pola magnetycznego na konwekcję naturalną płynu paramagnetycznego w geometrii sześciennej z różną konfiguracją ścianek chłodzonych i ogrzewanych. Badania nad konwekcją termomagnetyczną były wykonywane również w geometrii cylindrycznej najpierw przez Filara [50], [51], [52] a następnie przez Wróbla [53]. Fornalik [1] przedstawiła kompleksowe opracowanie zawierające termodynamiczny opis zależności występujących w przypadku płynu paramagnetycznego umieszczonego w niejednorodnym polu magnetycznym.. 1.5 Cel pracy. Na podstawie analizy danych literaturowych można stwierdzić, że nie została podjęta próba analizy wpływu silnego pola magnetycznego na obszar przejścia pomiędzy konwekcją laminarną a turbulentną, ponieważ wszystkie dotychczasowe badania koncentrowały się na obszarze konwekcji laminarnej. W. przedstawionej. pracy. konwekcja. termomagnetyczna. płynu. paramagnetycznego zostanie zbadana w geometrii sześciennej, w konfiguracji Rayleigha-Bénarda. Celem pracy jest zaproponowanie modelu matematycznego zjawiska konwekcji termomagnetycznej. Na podstawie zaproponowanego modelu matematycznego przeprowadzone zostanie. szereg. trójwymiarowych,. niestacjonarnych. symulacji. numerycznych. obejmujących obszar przejściowy i turbulentny konwekcji termomagnetycznej. Z uzyskanych wyników symulacji numerycznych zbudowana zostanie „mapa stabilności” wskazująca typ konwekcji termomagnetycznej na podstawie liczb podobieństwa. W celu weryfikacji zaproponowanego modelu matematycznego konwekcji termomagnetycznej oraz otrzymanych wyników numerycznych wykonana zostanie. 13.

(24) analiza eksperymentalna. Będzie miała na celu potwierdzenie obszarów konwekcji termomagnetycznej zidentyfikowanych na drodze analizy numerycznej. Zaprezentowane. zostaną. otrzymane. wyniki. eksperymentalne. wpływu. niejednorodnego pola magnetycznego na konwekcję termomagnetyczną w obszarze przejścia pomiędzy konwekcją laminarną a turbulentną. Dokonane zostanie porównanie przejścia pomiędzy obszarem przejściowym a turbulentnym dla konwekcji naturalnej (na podstawie danych literaturowych) i termomagnetycznej (na podstawie wyników eksperymentalnych). Przedstawiona zostanie korelacja liczby Nusselta i termomagnetycznej liczby Rayleigha.. 14.

(25) 2 Model matematyczny konwekcji termomagnetycznej. 2 Model matematyczny konwekcji termomagnetycznej 2.1 Magnetyzacja materiałów paramagnetycznych Equation Section 2 Zgodnie z prawem Curie [38] namagnesowanie próbki paramagnetycznej można mo wyrazić jako:. M =C. b. θ. (2.1). gdzie: θ to temperatura, natomiast C jest stałą materiałową zwanąą stałąą Curie. Implikacją powyższego ższego równania jest spadek wartości współczynnika współczynnik masowej podatności ci magnetycznej wraz ze wzrostem temperatury (równania (1.1))-(1.4)). W niskich temperaturach lub dla bardzo silnych silnych pól magnetycznych, namagnesowanie traci liniową liniow zależność od pola zewnętrznego i wykazuje nasycenie ( M max ) dla całkowitego uporządkowania dkowania dipoli w danej objętości. obj Zależność namagnesowania paramagnetyków od zewnętrznego pola magnetycznego, netycznego, znormalizowana do temperatury pomiaru wraz z prawem Curie i zakresem jego stosowalności stosowalno przedstawiono na rysunku 2−1.. Na rysunku 2−2 przedstawiono analogicznąą zależność zale obliczoną dla jonów Gd3+ występujących w stosowanych roztworach eksperymentalnych. eksperyment. Rys. 2-1 Zależność namagnesowania paramagnetyków od zewn zewnętrznego trznego pola magnetycznego, znormalizowana do temperatury pomiaru wraz z prawem Curie i zakresem jego stosowalności stosowalno w oparciu o dane literaturowe [3]. 15.

(26) Zakres stosowalności prawa Curie dla swobodnych jonów występujących w badanych roztworach Gd(NO3)3·6H2O w zależności od temperatury oraz zewnętrznego pola magnetycznego został również przedstawiony na rysunku 2−3. Z analizy (rys. 2−3) wynika, iż prawo Curie z powodzeniem opisuje namagnesowanie jonów Gd3+ w temperaturach powyżej θ > 100 [K] i polach magnetycznych mniejszych niż b 0 < 20 [T] .. Rys. 2-2 Zależność namagnesowania paramagnetyków od zewnętrznego pola magnetycznego, znormalizowana do temperatury dla jonów Gd3+[54]. Rys. 2-3 Zależność namagnesowania od przyłożonego zewnętrznego pola magnetycznego dla różnych temperatur dla swobodnych jonów Gd3+ [54]. 16.

(27) 2 Model matematyczny konwekcji termomagnetycznej. 2.2 Siły masowe Siła grawitacyjna. Konwekcja naturalna jest ruchem płynu wynikającym z oddziaływania siły grawitacyjnej na cząsteczki płynu posiadające pewną gęstość. Siła oddziaływania grawitacji na płyn może być zapisana zależnością: f gr = gρ .. (2.2). Jeżeli gęstość płynu zależy od temperatury to na cząsteczki posiadające różną temperaturę siła grawitacyjna działa różną siłą, czego efektem jest powstanie grawitacyjnej siły wyporu: f B = g( ρ − ρ0 ),. (2.3). gdzie: ρ 0 to gęstość w temperaturze referencyjnej. Model Boussinesqa przyjmuje, że gęstość jest wielkością stała, poza równaniem na siłę wyporu oraz jest on uzasadniony jedynie dla małych różnic temperatur [55]. Dzięki niemu, można wyrazić zależność gęstości płynu od temperatury, jako:. ρ ≃ ρ 0 1 − β (θ − θ 0 )  ,. (2.4). gdzie: θ0 to temperatura referencyjna, natomiast współczynnik rozszerzalności termicznej β zdefiniowany jest równaniem:. β =−. 1  ∂ρ ρ  ∂θ.  . . (2.5). Wprowadzając zależność (2.4) do równania opisującego grawitacyjną siłę wyporu działającą na jednostkę objętości płynu (2.3), można uzyskać wzór:. f B = −gρ0 β (θ − θ0 ) .. (2.6). Efektem działania grawitacyjnej siły wyporu, jest przemieszczanie się płynu zimnego zgodnie z kierunkiem i zwrotem siły grawitacyjnej. Ruch płynu ciepłego następuje w tym samym kierunku, ale ze zwrotem przeciwnym.. 17.

(28) Siła magnetyczna. Rozwój materiałów nadprzewodzących, których krytyczna temperatura była powyżej temperatury wrzenia ciekłego azotu, umożliwił budowę magnesów nadprzewodzących zdolnych generować silne pola magnetyczne (o indukcji magnetycznej do 10-20 [T]). Gradient tak silnego pola magnetycznego jest w stanie wpływać na przepływ płynów diamagnetycznych i paramagnetycznych, nieprzewodzących prądu elektrycznego. Siła działająca na paramagnetyk w polu magnetycznym na jednostkę objętości może być opisana następująco [56]: f mg = µ 0 ( M ⋅ ∇ ) H. (2.7). łącząc powyższy wzór z równaniem (1.1) otrzymuje się:. f mg = µ0 χ ( H ⋅∇ ) H =. µ0 χ 2. ∇H 2. (2.8). Wykorzystując relację (1.3), (1.4) oraz dlatego, że χ ≪ 1 , wzór może przyjąć postać: f mg ≈. χ ∇b 2 2 µ0. (2.9). Tagawa [43] wykorzystując: powyższe równanie, prawo Curie (równanie (2.1)) oraz aproksymację Boussinesqa (2.4) zaproponował równanie opisujące magnetyczną siłę wyporu:.  1  χ m 0 ρ0 β (θ − θ0 ) 2 f M = − 1 + ∇b  2 µ0  βθ0 . (2.10). gdzie: χ m 0 - masowa podatność magnetyczna w temperaturze referencyjnej θ 0 . Z równania (2.10) wynika, że płyn paramagnetyczny, którego temperatura θ jest niższa od temperatury referencyjnej θ0 jest przyciągany w kierunku kwadratu gradientu pola magnetycznego ∇b 2 . W przeciwnym przypadku, gdy θ > θ 0 płyn jest odpychany w przeciwną stronę.. 18.

(29) 2 Model matematyczny konwekcji termomagnetycznej. 2.3 Równania warunków. Równania opisujące proces konwekcji termomagnetycznej wynikają z zasady zachowania: pędu, masy i energii. Równanie pędu składa się z równania Naviera-Stokesa wraz z dodatkowymi członami opisującymi siłę wyporu grawitacyjnego f B oraz siłę wyporu magnetycznego f M i można przedstawić następująco:. ρ. ∂u + ( u ⋅ ∇ ) uρ = −∇ p + µ∇ 2u + f B + f M , ∂t. (2.11). f B = −gρ0 β (θ − θ0 ) ,. (2.12).  1  χ m 0 ρ0 β (θ − θ0 ) 2 f M = − 1 + ∇ b0 ,  2 µ0  βθ 0 . (2.13). gdzie: b 0 indukcja magnetyczna w centrum magnesu nadprzewodzącego, p - ciśnienie,. u - wektor prędkości, t - czas. Równanie energii:. ∂θ + ( u ⋅∇ ) θ = α∇2θ . ∂t. (2.14). ∂u = 0.. (2.15). Równanie ciągłości:. Do obliczenia rozkładu pola magnetycznego w rozpatrywanej przestrzeni wewnątrz cewki magnesu nadprzewodzącego wykorzystywane jest prawo Biota-Savarta:. b=. µ0 I ds × r ∫ 4π cewki r 3. (2.16). gdzie: I - natężenie prądu, ds - wektor elementu przewodnika z prądem, r - wektor wodzący, r - długość wektora wodzącego.. 19.

(30) 2.4 Liczby podobieństwa. Bezwymiarowe. liczby. wykorzystane. podobieństwa. do. opisu. konwekcji. termomagnetycznej mają następującą postać: •. termiczna liczba Rayleigha:. Ra T = •. gβ (θ h − θc ) D3. αν. ,. magnetyczna liczba Rayleigha [42]:  1   γ gβ (θ h − θ c ) D 3  Ra M =  1 +  , 2αν   βθ 0  . •. (2.17). (2.18). gamma - liczba opisująca magnetyzację płynu [42]:. χ b 0 max , γ= µ m gD 2. •. termomagnetyczna liczba Rayleigha: Ra TM =Ra T +Ra M ,. •. (2.20). liczba Prandtla:. Pr =. 20. (2.19). ν . α. (2.21).

(31) 3 Analiza numeryczna. 3 Analiza numeryczna 3.1 Modelowanie numeryczne Equation Section (Next) W celu rozwiązania równań: ciągłości, Naviera-Stokesa z dodatkowymi siłami wyporu oraz energii (równania (2.11)-(2.15)), zaprezentowanych w rozdziale 2.3, poddano je dyskretyzacji przy wykorzystaniu metody objętości skończonych – algorytmu SIMPLE (Semi-Implicit Method for Pressure-Linked Equations), powszechnie stosowanego w zagadnieniach związanych z przepływem płynu i wymianą ciepła [57][58]. Wykorzystano niejednorodną, zagęszczoną przy brzegach, przesuniętą siatkę, w której wielkości skalarne zdefiniowane zostały wewnątrz objętości skończonych, natomiast składowe prędkości zostały zdefiniowane na odpowiadających im ściankach bocznych. Do obliczenia niezbędnej gęstości siatki oraz rozmiaru kroku czasowego do symulacji przejścia pomiędzy konwekcją laminarną a turbulentną posłużono się modelem Kołmogorowa [59]. Zakładając, iż właściwości przepływu w małych skalach zależą wyłącznie od lepkości płynu i współczynnika dyssypacji można obliczyć charakterystyczne skale dla najmniejszych struktur - skale długości: 1/4. ν 3  η =  , ε . (3.1). ν  τη =   , ε . (3.2). oraz skalą czasu: 1/2. gdzie współczynnik dyssypacji lepkiej ε można oszacować, jako:. ε≈. u3 . D. (3.3). Po obliczeniu wymaganego rozmiaru siatki, porównano siatki 82x82x82 oraz 122x122x122 zagęszczone przy brzegu w celu zbadania wpływu rozmiaru siatki obliczeniowej na otrzymywane rezultaty. Dla wybranego przykładu - liczby Prandtla równej. 21.

(32) 584 różnica liczby Nusselta pomiędzy dwoma wybranymi siatkami wyniosła poniżej 1%. Ze względu na krótszy czas obliczeń wybrano siatkę o rozmiarze 82x82x82. Krok czasowy ustalono na poziomie ∆ t = 0,1 [s] . Obliczenia prowadzono dla całkowitego czasu obliczeń t > 3000 [s] oraz do momentu uzyskania założonej zbieżności dla równań: pędu, energii i. ciągłości, ustalonej na poziomie 10 −8 . Na rysunku 3−1a została przedstawiona siatka obliczeniowa składająca się z 82 objętości skończonych w każdym kierunku, wraz z zaznaczonym przekrojem centralnym, który został zaprezentowany na rysunku 3−1b.. Rys. 3-1 Siatka obliczeniowa: a) w całości wraz z zaznaczoną płaszczyzną centralną, b) tylko płaszczyzna centralna siatki. 3.2 Schematy numeryczne. Dobór schematów numerycznych jest istotnym elementem obliczeń, ponieważ wybór niewłaściwego może generować błędy numeryczne lub nierealistyczne rozwiązania. Na dokładność uzyskanych wyników wpływa również wykorzystanie schematów numerycznych wyższego rzędu. Zastosowanie schematów drugiego rzędu jest uważane za zupełnie wystarczające dla zastosowań inżynieryjnych [60]. W celu zachowania dokładności. 22.

(33) 3 Analiza numeryczna rozwiązania, należy dyskretyzować człon konwekcyjny, dyfuzyjny, źródła i człon niestacjonarny wykorzystując ten sam rząd stosowanych schematów numerycznych. W przeciwnym razie rozwiązanie numeryczne będzie cechować się dokładnością najmniej dokładnego członu równania [61]. W przypadku symulacji obszaru przejściowego pomiędzy konwekcją laminarną, a turbulentną, niepożądanym schematem jest np. schemat UDS (Upwind-Difference Scheme)[62], charakteryzujący się dużą stabilnością i szybką zbieżnością kosztem sporej dyfuzji numerycznej (modelowane medium cechuje się wyższą dyfuzją niż w rzeczywistości) oraz zależności od kierunku przepływu płynu (pochodna aproksymowana jest ilorazem różnicowym „wstecz” albo „w przód” w zależności od kierunku przepływu). Dlatego, w rozważanym zagadnieniu, wykorzystany został schemat CDS (Central-Difference Scheme) drugiego rzędu [63], pozbawiony powyższych wad. Uwzględniając przedstawione uwagi, do rozwiązania członu niestacjonarnego, wykorzystany został niejawny schemat aproksymacji drugiego rzędu wykorzystujący trzy kolejne kroki czasowe (Implicit Three-Consecutive Time-Steps Integration Scheme)[63].. 3.3 Właściwości płynu. Do symulacji numerycznych wykorzystano dane literaturowe [64] własności płynu paramagnetycznego odpowiadające 80% masowemu wodnemu roztworowi gliceryny z dodatkiem 0,8 mol Gd ( NO3 )3 ×6H2O , którego właściwości fizyczne zawarto w tabeli 3-1. Własności płynu oznaczone gwiazdką Bednarz [64] oszacował na podstawie danych literaturowych. W celu otrzymania różnej liczby Prandtla w obliczeniach skalowano jedynie dynamiczny współczynnik lepkości, pozostałe własności utrzymując niezmienione. Założono, zgodnie z wcześniejszą analizą [65], iż taki płyn można uznać za: •. newtonowski,. •. nieściśliwy,. •. paramagnetyczny,. •. nieprzewodzący prądu elektrycznego.. 23.

(34) Tabela 3-1 Właściwości fizyczne płynu w temperaturze 298 [K], zastosowanego do obliczeń numerycznych. Własność. Symbol. Wartość. Dyfuzyjność termiczna*. α [m 2 /s]. 1, 01 ⋅10 −7. Gęstość. ρ [kg/m 3 ]. 1463. Dynamiczny współczynnik lepkości. µ [kg/m ⋅ s]. 8, 69 ⋅10 −2. Kinematyczny współczynnik lepkości. ν [m 2 /s]. 5,9 ⋅10−5. Współczynnik przewodzenia ciepła*. λ [W /m ⋅ K ]. 0, 397. Współczynnik rozszerzalności termicznej. β [1 / K]. 5, 2 ⋅10 −4. Masowa podatność magnetyczna. χ m [m 3 /kg]. 2, 309 ⋅10 −7. Liczba Prandtla. Pr [-]. 584. * właściwości płynu oszacowane przez Bednarza na podstawie literatury [64].. 3.4 Rozkład pola magnetycznego. Rozkład pola magnetycznego i jego gradient zostały obliczone dla rzeczywistego elektromagnesu, na którym przeprowadzono eksperymenty – Superconducting Helium Free Magnet HF10-100VHT-B Sumitomo Heavy Industries. Jest to magnes składający się z dwóch cewek: z wewnętrznej, wykonanej z Nb3Sn oraz zewnętrznej, z NbTi. Ich schemat został przedstawiony na rysunku 3−2. Rozkład pola magnetycznego został obliczony poprzez scałkowanie równania Biota-Savarta (równanie (2.16)) dla rzeczywistych wymiarów, każdej z cewek magnesu. Na rysunku 3−3 przedstawiono wyniki obliczeń rozkładu pola magnetycznego oraz jego gradientu. Dodatkowo na rysunku 3−3 została oznaczona pozycja układu eksperymentalnego wykorzystana w obliczeniach numerycznych. Na rysunku 3−4 przedstawiono kontury natężenia gradientu kwadratu indukcji magnetycznej badanego sześcianu przy | b 0 |max =10 [T]. Gradient w rozpatrywanym układzie osiągał wartość do 900 [T2/m]. Zmianę rozkładu pola oraz gradientu kwadratu indukcji magnetycznej w centrum sześcianu ( yc =0,016 [m], xc =0,016 [m], 0 ≤ zc ≤ 0,032 [m]) w całym zakresie stosowanej indukcji magnetycznej ( 0 [T] ≤ b 0 rysunku 3−5.. 24. max. ≤ 10 [T] ) przedstawiono na.

(35) 3 Analiza numeryczna. Rys. 3-2 Schemat przedstawiający przedstawiają przekrój poprzeczny przez cewki magnesu HF10-100VHT-B HF10. Rys. 3-3 Pozycja zycja układu eksperymentalnego naniesiona na: a) rozkład rozkład pola magnetycznego dla |b0|max=10 [T], b) gradient kwadratu indukcji magnetycznej dla |b0|max=10 [T] [66]. 25.

(36) Rys. 3-4 Kontur natężenia gradientu kwadratu indukcji magnetycznej badanego sześcianu przy |b0|max=10 [T] [67]. Rys. 3-5 Zmiana: a) rozkładu pola magnetycznego oraz b) gradientu kwadratu indukcji magnetycznej we wnętrzu sześcianu (xc=0,016 [m], yc=0,016 [m], 0 [m]≤zc≤0,032 [m]) dla całego rozpatrywanego zakresu indukcji magnetycznej (0 [T]<|b0|max<10 [T]) [67]. 26.

(37) 3 Analiza numeryczna. 3.5 Geometria i warunki brzegowe. Do obliczeń numerycznych, tak samo jak do eksperymentu została wybrana geometria du kartezjańskiego kartezja D = 0, 032 [m] . Środek układu. sześcienna. Długość ść boku wynosiła. umiejscowiony został na rogu dolnej ścianki sześcianu cianu jak przedstawiono na rysunku 3−6. Środek sześcianu położony żony został. z = 0,11 [m]. ponad centrum zewnętrznego zewn pola. magnetycznego | b 0 |max . Temperaturę Temperatur ogrzewanej ścianki cianki dolnej ustalono na θ h = 25 [ o C] , natomiast chłodzonej ścianki górnej na θ c = 20 [ o C] .. Rys. 3-6 Geometria układu ukł wraz z warunkami brzegowymi. Warunki brzegowe modelowanego układu: •. Dla wszystkich ścianek geometrii sześciennej: sze u = 0,. •. Ścianki boczne sąą adiabatyczne, adiabatyczne odpowiednio:. •. Horyzontalna ścianka dla zc = 0 [m] : θ = θ h ,. •. Horyzontalna ścianka dla zc = D = 0, 032 [m] : θ = θ c .. ∂θ ∂θ = 0, = 0, ∂Y ∂X. Warunki początkowe tkowe dla obliczeń oblicze numerycznych: •. Dla t = 0 : u = 0, θ = (θ h − θc ) / 2. 27.

(38) Wymianę ciepła opisano przy pomocy liczby Nusselta, jako stosunek ciepła wymienianego na drodze konwekcji Q conv , do strumienia ciepła wymienianego na drodze przewodzenia Qcond : Nu=. Qconv Qcond. (3.4). W obliczeniach numerycznych, liczbę Nusselta obliczono na ściance ciepłej ( zc = 0 ) zgodnie z równaniem:  ∂θ convection   dXdY ∂Z  Z =0 0 0 = DD  ∂θ   dXdY ∫0 ∫0  conduction ∂Z Z =0 DD. Nu hot. ∫ ∫ . (3.5). 3.6 Wprowadzenie do wyników obliczeń Symulacje numeryczne przeprowadzono dla liczb Prandtla w zakresie 5, 8 < Pr < 584. Poniżej przedstawiono wyniki dla przypadków Pr=20, 100, 146 i 584 w zakresie 0 [T] ≤ b 0 b0. max. max. ≤ 10 [T] . Obliczenia numeryczne wykonano jedynie dla całkowitych wartości. przy kroku równym 1 [T]. W przedstawionych wynikach wykorzystane jest pojęcie indukcji magnetycznej w. centrum magnesu nadprzewodzącego b 0. max. , zamiast parametru gamma (γ ) opisującego. magnetyzację płynu, w celu wygodniejszego porównania wyników numerycznych z wynikami uzyskanymi na drodze eksperymentu. Wartości parametru gamma odpowiadające b0. max. dla symulowanego płynu posiadającego podatność χ m = 2, 309 ⋅ 10 −7 zebrano w tabeli. 3-2. Tabela 3-2 Wartości indukcji magnetycznej w centrum magnesu nadprzewodzącego i odpowiadające mu wartości parametru gamma dla symulowanego płynu. b0. γ. 28. max. 0. 1. 2. 3. 4. 0,00. 0,59. 2,34. 5,27. 9,37. 5. 6. 7. 8. 9. 10. 14,64 21,08 28,70 37,48 47,44 58,56.

(39) 3 Analiza numeryczna. 3.7 Wariant obliczeń nr 1 (Pr=20) Obliczenia zostały przeprowadzone dla liczby Prandtla: Pr=20, przy ∆ θ = 5 [K] , co razem z własnościami płynu odpowiadało termicznej liczbie Rayleigha: Ra T = 4, 06 ⋅106 . W przypadku. | b 0 |max = 0 [T]. (rys. 3−7), można zidentyfikować szereg komórek. konwekcyjnych, które przy | b 0 |max = 1 [T] (rys. 3−8) zostają znacząco zniekształcone. W całym pozostałym zakresie zastosowanej indukcji magnetycznej 2 [T] ≤| b 0 |max ≤ 10 [T] brak jest widocznej organizacji pól temperatury i lini prądu (rys. 3−9 – rys. 3−17), a ich chaotyczność postępuje wraz ze wzrostem indukcji magnetycznej. Zmniejsza się warstwa przyścienna, a centralna część płynu osiąga niemal jednorodną temperaturę. Obserwuje się wtrącenia zarówno z warstwy górnej chłodzonej, jak i dolnej ogrzewanej, które wraz ze wzrostem indukcji magnetycznej w centrum magnesu sięgają coraz dalej, osiągając przy tym większą prędkość. Jednak, przy | b 0 |max = 10 [T] wtrącenia stają się bardzo krótkie i szybkie. W tym przypadku płyn posiada niemal jednorodną temperaturę.. 29.

(40) a). b). c). d) Z. X. Z. Y. X. Y. e) Y. Z. X. Rys. 3-7 Wyniki symulacji dla RaTM=4,06·106, Pr=20, przy |b0|max=0 [T]: a) izopowierzchnie temperatury, b) linie prądu wraz z polem prędkości w przekroju centralnym x-y, c-e) pole temperatury – przekrój centralny w osi: z-x, z-y i y-x. 30.

(41) 3 Analiza numeryczna a). b). c). d) Z. Y. Z. X. X. Y. e) Y. Z. X. Rys. 3-8 Wyniki symulacji dla RaTM=1,29·107, Pr=20, przy |b0|max=1 [T]: a) izopowierzchnie temperatury, b) linie prądu wraz z polem prędkości w przekroju centralnym x-y, c-e) pole temperatury – przekrój centralny w osi: z-x, z-y i y-x. 31.

(42) a). b). c). d) Z. Y. Z. X. X. Y. e) Y. Z. X. Rys. 3-9 Wyniki symulacji dla RaTM=3,95·107, Pr=20, przy |b0|max=2 [T]: a) izopowierzchnie temperatury, b) linie prądu wraz z polem prędkości w przekroju centralnym x-y, c-e) pole temperatury – przekrój centralny w osi: z-x, z-y i y-x. 32.









(51) 3 Analiza numeryczna. 3.8 Wariant obliczeń nr 2 (Pr=100) Obliczenia zostały przeprowadzone dla liczby Prandtla: Pr=100, przy ∆ θ = 5 [K] , co razem z własnościami płynu odpowiadało termicznej liczbie Rayleigha: Ra T = 8,15 ⋅105 . W przypadku | b 0 |max = 0 [T] , można zaobserwować brak symetrii komórek konwekcyjnych na wizualizacjach izopowierzchni temperatury (rys. 3−18a), przekrojach pól temperatury (rys. 3−18cde) oraz liniach prądu (rys. 3−18b). Po przyłożeniu słabego pola magnetycznego 1 [T] ≤| b 0 |max ≤ 3 [T]. następuje. uporządkowanie przepływu płynu oraz uporządkowanie pól temperatury, które stają się w pełni symetryczne względem centralnej pionowej osi przekroju (rys. 3−19 – rys. 3−21). Podczas dalszego zwiększania indukcji magnetycznej, przy | b 0 |max ≥ 4 [T] można zaobserwować narastające nieuporządkowanie w polach temperatury (rys. 3−22 – rys. 3−28acde) na skutek intensywniejszego mieszania się płynu poprzez wtrącenia zarówno z warstwy górnej chłodzonej, jak i dolnej ogrzewanej oraz przemieszczanie się głównego rdzenia przepływu względem osi centralnej geometrii sześciennej (rys. 3−22 – rys. 3−28b).. 41.












(63) 3 Analiza numeryczna. 3.9 Wariant obliczeń nr 3 (Pr=146) Obliczenia zostały przeprowadzone dla liczby Prandtla: Pr=146, przy ∆ θ = 5 [K] , co razem z własnościami płynu odpowiadało termicznej liczbie Rayleigha: Ra T = 5, 61 ⋅105 . W przypadku braku zewnętrznego pola magnetycznego ( | b 0 |max = 0 [T] ) oraz przyłożonego słabego pola magnetycznego ( | b 0 |max = 1 [T] ) można zaobserwować brak symetrii względem pionowej osi centralnej sześcianu, zarówno w przypadku pól temperatury (rys. 3−29 – rys. 3−30cde), jak i liniach prądu (rys. 3−29 – rys. 3−30b). Po przyłożeniu słabego pola magnetycznego 2 [T] ≤| b 0 |max ≤ 4 [T] następuje uporządkowanie przepływu płynu (rys. 3−31 – rys. 3−33b) oraz uporządkowanie pól temperatury (rys. 3−31 – rys. 3−33cde), które stają się w pełni symetryczne. Podczas dalszego zwiększania indukcji magnetycznej, przy | b 0 |max ≥ 5 [T] można zaobserwować narastające nieuporządkowanie w polach temperatury. (rys. 3−34 – rys. 3−39cde) na skutek intensywnego mieszania się płynu poprzez wtrącenia zarówno z warstwy górnej chłodzonej, jak i dolnej ogrzewanej oraz przemieszczanie się głównego rdzenia przepływu względem osi centralnej geometrii sześciennej (rys. 3−34 – rys. 3−39b).. 53.

(64) a). b). c). d) Z. Y. Z. X. X. e) Y. Z. X. Rys. 3-29 Wyniki symulacji dla RaTM=5,61·105, Pr=146, przy |b0|max=0 [T]: a) izopowierzchnie temperatury, b) linie prądu wraz z polem prędkości w przekroju centralnym x-y, c-e) pole temperatury – przekrój centralny w osi: z-x, z-y i y-x. 54. Y.




(68) a). b). c). d) Z. Y. Z. X. X. e) Y. Z. X. Rys. 3-33 Wyniki symulacji dla RaTM=2,01·107, Pr=146, przy |b0|max=4 [T]: a) izopowierzchnie temperatury, b) linie prądu wraz z polem prędkości w przekroju centralnym x-y, c-e) pole temperatury – przekrój centralny w osi: z-x, z-y i y-x. 58. Y.







(75) 3 Analiza numeryczna. 3.10 Wariant obliczeń nr 4 (Pr=584). Obliczenia zostały przeprowadzone dla liczby Prandtla: Pr=584, przy ∆ θ = 5 [K] , co razem z własnościami płynu odpowiadało termicznej liczbie Rayleigha: Ra T = 1, 40 ⋅105 . W przypadku | b 0 |max = 0 [T] , można zaobserwować brak symetrii zarówno w przypadku pól temperatury (rys. 3−40cde), jak i linii prądu (rys. 3−40b). Po przyłożeniu słabego pola magnetycznego 1 [T] ≤| b 0 |max ≤ 2 [T] następuje uporządkowanie przepływu płynu (rys. 3−41 – rys. 3−42b) oraz uporządkowanie pól temperatury (rys. 3−41 – rys. 3−42cde), które stają się w pełni symetryczne. Podczas dalszego zwiększania indukcji magnetycznej, przy | b 0 |max = 3 [T] wkrada się asymetria w polach temperatury (rys. 3−43cde) oraz przepływie. (rys. 3−43b), jednak przy | b 0 |max ≥ 4 [T] całość nadal wykazuje silną symetrię (rys. 3−44 – 3−49). Przy | b 0 |max = 10 [T] można zaobserwować pewne zaburzenia w polach temperatury (rys. 3−50cde), spowodowane drobnymi wtrąceniami zarówno z warstwy górnej chłodzonej, jak i dolnej ogrzewanej.. 65.












(87) 3 Analiza numeryczna. 3.11 Analiza porównawcza wyników numerycznych. Na rysunku 3−51 zebrano termiczne wyniki symulacji numerycznej przedstawione w postaci liczby Nusselta, obliczonej ścianie ciepłej ( zc = 0 , zgodnie z równaniem (3.5)), dla symulowanego płynu dla wybranych liczb Prandtla wynoszących: 20, 100, 146 i 584. Wyniki przedstawiono, jako różnicę liczby Nusselta dla danej indukcji magnetycznej w centrum magnesu ( | b 0 |max ) Nu b i liczby Nusselta w przypadku braku pola magnetycznego Nu 0 . Intensyfikacja wymiany ciepła dzięki wykorzystaniu konwekcji termomagnetycznej, w odniesieniu do konwekcji naturalnej [ (Nu b /Nu 0 ) − 1 ] wacha się od 15% ( | b 0 |max = 1 [T] ) do 330% ( | b 0 |max = 10 [T] ).. Rys. 3-51 Termiczne wyniki symulacji numerycznej, przedstawione jako: różnica liczby Nusselta dla danej indukcji magnetycznej w centrum magnesu i liczby Nusselta w przypadku braku pola magnetycznego dla różnych wartości indukcji magnetycznej w centrum magnesu. Na rysunku 3−52 przedstawiono liczbę Nusselta w odniesieniu do termomagnetycznej liczby Rayleigha. Tło rysunku określa typ konwekcji: kolor szary - obszar stabilny, biały obszar turbulentny. W wynikach numerycznych zidentyfikowano destabilizację wyników: 77.

(88) pól temperatury oraz prądów płynu, które dla płynów o liczbie Prandtla równej: 20, 100, 146 i 584, pojawiły się odpowiednio przy: | b 0 |max = 2 [T] , | b 0 |max = 4 [T] , | b 0 |max = 5 [T] i | b 0 |max = 10 [T] . Na podstawie wartości termomagnetycznej liczby Rayleigha w tych. przypadkach, można przyjąć, że to punkt przejścia pomiędzy obszarem stabilnym i turbulentnym w obliczeniach numerycznych wynosi Ra TM = 3 ⋅107 .. Rys. 3-52 Wyniki symulacji numerycznej - liczba Nusselta w odniesieniu do termomagnetycznej liczby Rayleigha. Heslot [9], przeprowadził analizę statystyczną dla konwekcji naturalnej, którą podzielił na podstawie termicznej liczby Rayleigha, na trzy obszary: przejściowy dla Ra T < 105 , lekkiej turbulencji dla 105 ≤ Ra T ≤ 4 ⋅107 oraz silnej turbulencji dla Ra T > 4 ⋅107 . Porównując otrzymane wyniki numeryczne konwekcji termomagnetycznej z analizą Heslota [9], można zauważyć, że:. •. wszystkie przeprowadzone symulacje numeryczne, w przypadku braku indukcji magnetycznej w centrum magnesu ( | b 0 |max = 0 [T] ) znajdują się w obszarze lekkiej turbulencji,. 78.