MODELOWANIE OSCYLACJI TEMPERATURY POWIERZCHNI GRZEJNEJ
WE WRZENIU W MIKROKANALE
Hubert Grzybowski
1a, Romuald Mosdorf
1b1Katedra Mechaniki i Informatyki Stosowanej, Wydział Mechaniczny, Politechnika Białostocka e-mail: agrzybowskihubert@wp.pl, bmosdorf@ii.pb.bialystok.pl
Streszczenie
W pracy analizowano wymianę ciepła i masy w układzie zawierającym mikrokanał. Przedstawiono model pozwala- jący symulować oscylacje wielkości opisujących wymianę ciepła i masy w mikrokanale (temperatury powierzchni grzejnej, ciśnienia, natężenia przepływu cieczy). Do opisu przepływu mieszaniny dwufazowej zastosowano jedno- wymiarowy model przepływu typu – ‘homogeneous flow model’. Model przepływu ciepła w ściance mikrokanału oparto na jednowymiarowym równaniu przewodnictwa ciepła z uwzględnieniem zmiennego w czasie strumienia ciepła pobieranego przez wrzącą ciecz. Otrzymane wyniki symulacji pokazują, że pojawienie się oscylacji parame- trów opisujących wymianę ciepła i masy w mikrokanale określone jest poprzez wzajemną interakcję pomiędzy układem zasilającym wymiennik ciepła w ciecz, a układem dostarczającym ciepło do mikrokanału. Uzyskany w modelu scenariusz pojawiania się oscylacji parametrów opisujących wymianę ciepła i masy obserwowany był w badaniach eksperymentalnych.
MODELLING OF MICROCHANNEL HEATING SURFACE TEMEPRATURE OSCILATIONS IN BOILING
Summary
In this paper the heat and mass transfer inside the systems with microchannel was analyzed. Proposed in the pa- per model allows to simulate the oscillations of parameters describing the heat and mass transfer inside microchannel (heating surface temperature, pressure, liquid flow rate). To describe the two-phase mixture flow the one-dimensional flow model type called 'homogeneous flow model' was used. The model of the heat transfer inside the microchannel walls was based on the one-dimensional heat conduction equation with changing in time the heat flux absorbed by the boiling liquid. The obtained simulation results show that the appearance of the oscilla- tions of parameters describing the heat and mass transfer inside the microchannel is determined by the mutual in- teractions between the liquid supply and the heat supply systems of the microchannel. The similar oscillations patterns of parameters describing the heat and mass transfer obtained from the model was observed in experi- mental researches.
1. WSTĘP
Wymienniki ciepła wyposażone w mikrokanały znalazły zastosowanie w mikroreaktorach chemicznych, mikro- mieszalnikach i mikrowymiennikach ciepła oraz w chłodzeniu procesorów. Mikrowymiennik potrafią odprowadzać duże ilości ciepła, ponieważ stosunek ich powierzchni grzejnej do objętości czynnika chłodzącego jest bardzo duży. Oscylacje natężenia przepływu, ciśnie- nia oraz temperatury mogą spowodować wibracje wy-
miennika, zaburzyć przepływ ciepła, co w szczególnych przypadkach mogą prowadzić do awarii systemu. Dlate- go zjawisko oscylacji temperatury powierzchni grzejnej, ciśnienia oraz natężenia przepływu zachodzące podczas dwufazowego przepływu cieczy przez mikrokanały jest zjawiskiem niepożądanym. Stało się ono przedmiotem badań wielu naukowców. Badania te omówione są między innymi w pracach: [1 ÷ 4]. Nieodpowiedni dobór
parametrów pracy systemu takich jak natężenie prze- pływu, ciśnienie, wymiary geometryczne, temperatura pracy może spowodować pojawienie się różnego rodzaju niestabilności. Niestabilności dzieli się na statyczne i dynamiczne. Niestabilności statyczne charakteryzuje się tym, że układ wytrącony z pierwotnego stanu rów- nowagi zmierza asymptotycznie do nowego stanu rów- nowagi. W przypadku niestabilności dynamicznych układ po wprowadzeniu zaburzenia próbuje osiągnąć nowy stan równowagi, lecz siły działające na układ powodują, że oscyluje on jedynie wokół stanu równowa- gi. Głównymi czynnikami odpowiedzialnymi za powsta- wanie niestabilności dynamicznych są opóźnienia czaso- we w propagacji zjawisk oraz sprzężenia zwrotne między zjawiskami zachodzącymi podczas przepływów dwufa- zowych. Wyróżnia się następujące typy niestabilności dynamicznych [4]:
1. oscylacje związane ze zmianą gęstości mieszaniny dwufazowej (density-wave oscillations),
2. oscylacje ciśnienia (pressure-drop oscillations), 3. oscylacje akustyczne (acoustic oscillations), 4. oscylacje termiczne (thermal oscillations).
Przeprowadzono wiele doświadczeń, w których zaobser- wowano niestabilności dynamiczne [6]. Wu i Cheng zaobserwowali w ogrzewanym mikrokanale oscylacje temperatury oraz ciśnienia o wysokiej amplitudzie [7 ÷ 9]. W badaniach niestabilności wrzenia w mikroka- nałach opisanych w pracy [5] obserwowano dynamiczną niestabilność pracy wymiennika. Mechanizm oscylacji polegał na cyklicznych zmianach zarówno temperatury powierzchni grzejnej jak i ciśnienia cieczy. Zaobserwo- wano, że powstająca w wyniku wrzenia w mikrokanale para blokuje przepływ cieczy, co prowadzi do wzrostu ciśnienia cieczy na wlocie do mikrokanału. W wyniku tego wzrostu następuje zwiększenie się natężenia prze- pływu cieczy. Proces ten prowadzi do zwiększenia ilości ciepła odbieranego przez wrzącą ciecz. W rezultacie następuje spadek temperatury powierzchni grzejnej i wrzenie staje się mniej intensywne lub zanika całkowi- cie. Pęcherze przestają się tworzyć i blokować przepływ, cykl powtarza się od nowa. Podobny mechanizm oscyla- cji parametrów opisujących wymianę ciepła i masy w mikrokanale obserwowano w badaniach eksperymen- talnych opisanych w pracy [4].
Najczęściej stosowanymi modelami przepływów dwufa- zowych są modele, w których przepływ dwóch faz - ciekłej i gazowej traktowany jest jak przepływ jednej
fazy o określonych parametrach (Homogeneous Flow Model) oraz modele, w których uwzględnia się stosunek ilości fazy ciekłej do gazowej w formie współczynnika modyfikującego właściwości mieszany dwufazowej (Separated Flow Model). Stosowany jest także model typu Drift-flux. Model ten wyróżnia się tym, że równa- nia zachowania masy, pędu oraz energii zapisywane są oddzielnie dla każdej z faz.
W pracy do opisu przepływu dwufazowego w mikroka- nale zastosowano model przepływu dwufazowego typu Homogeneous Flow Model. Zaproponowany model pozwala na symulację oscylacji temperatury powierzchni grzejnej, natężenia przepływu cieczy oraz spadku ciśnie- nia wzdłuż mikrokanału. W pracy przedstawiono wy- brane symulacje, które uzasadniają wyniki pomiarów przedstawionych w pracy [4].
2. MODEL WYMIANY CIEPŁA I MASY W MIKROKANALE
Rozpatrywano wymianę ciepła i masy zachodzącą we fragmencie mikrokanału pokazanym na rys. 1. Strzałką oznaczono kierunek przepływu mieszaniny dwufazowej.
Nad elementem i o wymiarach δ x δ przepływa wrząca ciecz wraz z natężeniem przepływu Gn. Do elementu i dostarczany jest strumień ciepła o natężeniu qd. Ciepło przekazywane jest przez ściankę kanału o grubości δ1 i odbierane przez płynącą przez kanał ciecz.
Rys. 1. Model układu, w którym analizowano wymianę ciepła i masy w mikrokanale
Równanie bilansu ciepła w elemencie i ma postać:
+ −
⋅
⋅
⋅
⋅
⋅
=
∆
−
⋅
−
∆
−
⋅ +
∆
−
⋅
−
∆
−
⋅ +
∆
⋅
⋅
⋅
−
∆
⋅
⋅
⋅
k Ti k Ti c F t
k T Ti k t
Ti TB
R t k T Ti k t
Ti TL k t
Ti k qu k t
qd
1 1
1 1
1 1
ρ δ δ δ δ
δ λδ δ
δ λδ
δ δ λδ δ
δ λδ δ δ δ
δ
(1) gdzie:
Ti jest średnią temperaturą elementu ‘i’,
( )
T tquk ik ⋅
δ
⋅δ
⋅∆ jest ciepłem odbieranym od elementu‘i’ w przedziale czasu
∆ t
,q
dk⋅ δ ⋅ δ ⋅ ∆ t
jest ciepłem dostarczanym do elementu ‘i’ w przedziale czasu∆ t
oraz
(
ik)
k
i T
T c⋅ −
⋅
⋅
⋅
⋅
δ δ
1ρ
+1δ
jest przyrostem energiielementu ‘i’.
Równanie (1) pozwala wyznaczyć temperaturę ścianki kanału, w kolejnych chwilach czasu (co przedział czasu
∆t), w postaci:
+ + + +
−
− +
+ =
TF TB TR TL k Ti k
Ti k qu k qd A k Ti k
Ti 4
1 1
1 1
λ δ δ
δ λ
δ δ
δ
(2) gdzie:
Fo t
A 2 c 1
∆ =
⋅
= ⋅ λ
ρ δ
Ciepło pobierane przez wrzącą ciecz zapisano w postaci:
) (Ti Tsat
qu=α⋅ − (3)
gdzie :
Ti– temperatura ścianki mikrokanału, Tsat– temperatura nasycenia,
α – współczynnik przejmowania ciepła.
Współczynnik przejmowania ciepła wyznaczono z kore- lacji Warriera [10], która została opracowana dla wrze- nia zachodzącego przy różnych temperaturach cieczy na wlocie (26, 40, 60 °C), przy przepływach cieczy w zakresie 557 ÷ 1600 kg/m2s i strumieniu ciepła do- starczanego do układu 0 ÷ 5.99 W/m2. Zgodnie z kore- lacją Warriera współczynnik α ma postać:
) ( 4
3 α
α E sp
Nu
= Nu (4)
gdzie:
d Nu k
h f
sp= 3
α ,
65 . 16 0
1
) 855 1 ( 3 . 5 6 0 .
1 Bo Bo x
E= + − − ,
), 0 . 2 361 . 5 814 . 5 767 . 3 883 . 1 1 ( 235 .
8 2 3 4 5
3= − β+ β − β + β − β
Nu
), 186 . 0 058 . 1 477 . 2 085 . 3 042 . 2 1 ( 235 .
8 2 3 4 5
4= − β+ β − β + β − β
Nu
Nu– liczby Nusselta,
αsp– współczynnik przejmowania ciepła dla przepływu jednofazowego,
Bo– liczba wrzenia, x – stopień suchości,
β – stosunek wysokości do szerokości kanału.
Przyjęto homogeniczny model przypływu, w którym przepływ dwufazowy traktuję się jako jednorodną mieszaninę fazy ciekłej i gazowej. W tym przypadku prędkości obydwu faz są równe. Przepływ mieszaniny dwufazowej w mikrokanale opisano z zastosowaniem jednowymiarowego równanie bilansu pędu dla przepły- wu wzdłuż kanału [5] w postaci:
) 0 ( 1 ) ( )
( 2
=
+
∂ + ∂
∂ +∂
∂
∂
Fvisc
GA z A z PA t
GA
ρ (5)
gdzie:
G – strumień masy, A– pole przekroju kanału,
Fvisc – siła tarcia spowodowana lepkością płynu.
z – odległość od początku kanału,
Dzieląc obie strony równania (5) przez A, po przekształ- ceniach otrzymano:
Fvisc
G z z P t
G −
∂
− ∂
∂
−∂
∂ =
∂
ρ
2
(6)
zakładając, że na długości kanału cieśnienie zmienia się liniowo równanie (6) zapisano w postaci:
L G P G L L
P P t
dG f
out in out
in ∆
−
−
− −
∂ = ρ ρ
2
1 2 (7)
gdzie:
∫
=
∆Pf 0LFviscdz – spadek ciśnienia wynikający z tarcia w kanale [11],
out
in P
P − – spadek ciśnienia nałożony na układ,
out in
G G
ρ ρ
2 2
− – spadek ciśnienia wynikający z przyspie-
szenia cieczy.
Zapisując różnicę ciśnień na początku i na końcu kanału w postaci Pin−Pout =∆P oraz rozpatrując zmiany wydatku w kanale w krótkim przedziale czasu, równanie (10) zapisano w postaci:
L G P
G L L
P t
G
G f
ou n in n n
n ∆
−
−
∆ +
∆ =
+ −
ρ ρ
2 2
1 1
(8)
Równanie (8) pozwala na wyznaczenie wartości Gn+1 w postaci:
n f
ou n in n
n G
L t G P
G L
t L
t
G P ∆ ⋅∆ +
−
∆ −
∆ +
⋅
=∆
+ ρ ρ
2 2
1 (9)
Spadek ciśnienia w kanale podczas przepływu dwufazo- wego opisany jest zależnością [12]:
g a ft
f P P P
P =∆ +∆ +∆
∆ (10)
gdzie:
Pft
∆ – spadek ciśnienia wywołany tarciem, Pa
∆ – spadek ciśnienia wywołany przyspieszeniem płynu,
Pg
∆ – spadek ciśnienia hydrostatycznego Zachodzi [6]:
h tp tp
f d
L G P f
ρ 2 2
=
∆ (11)
gdzie:
ρtp –gęstość mieszaniny dwufazowej:
tp v tp l
tp ρ ε ρ ε
ρ = (1− )+ [12]
ρl– gęstość cieczy ρv – gęstość pary εtp – stopień zapełnienia
+ −
=
l v tp
x x
ρ ε ρ
1 1
1 [6]
ftp– współczynnik tarcia, zachodzi [6]: 0.25 Re
079 .
= 0
ftp , zaś:
tp
Gdh
= µ Re
µtp – współczynnik lepkości mieszaniny dwufazowej:
l v
tp xµ xµ
µ = +(1− ) , [11]
Spadek ciśnienia spowodowany przyspieszeniem zapisa- no w postaci [11]:
dz d G
z
Pa P a tp
∂ =
= ∂
∆ ρ
)
( (12)
W dalszych rozważaniach pominięto wielkość ∆Pg. Strumień masy Gn+1, opisany zależnością (9), jest funk- cją współczynnika wypełnienia kanału parą xRównanie (4) po przekształceniach pozwala na wyznaczenie ter- modynamicznego stopnia suchości mieszaniny dwufazo- wej x w funkcji współczynnika przejmowania ciepła:
65 . 0
1
16 1
3 4
) 855 1 ( 3 . 5
6 1
⋅
−
⋅
−
⋅
−
⋅ −
⋅
= Bo
Nu Bo Nu x αsp
α
(13)
Równania (2) i (9) pozwalają na wyznaczenie zmian temperatury powierzchni grzejnej oraz masowego wy- datku cieczy w kolejnych chwilach czasu. Schemat obliczeń pokazano na rys. 2. Cykl rozpoczyna się od obliczenia temperatury ścianki mikrokanału Ti na podstawie równania (1). Zmiana temperatury ścianki grzejnej powoduję zmianę ciepła odbieranego qu, co powoduje zmienię stopnia zapełnienia kanału x, a w konsekwencji dochodzi do zmiany natężenia przepływu G. Nowy przepływ charakteryzuje się nową wartością współczynnika przejmowania ciepła α. Temperatura ścianki grzejnej Ti oraz nowa wartość współczynnika α określa nową wartość strumienia ciepła odbieranego qu. Cały cykl powtarza się od początku, gdyż nowa wartość qu zmienia temperaturę ścianki kanału Ti.
Rys. 2.Schemat obliczeń zastosowany w modelu.
Obliczenia wykonano dla wody przy średnicy hydrau- licznej mikrokanału dh = 100 μm, natężeniu przepływu G = 200 kg.m-2s-1 oraz stopniu suchości mieszaniny duwfazowej x = 0.6. Symulowano oscylacje parametrów opisujących badany układ dla różnych gęstości strumie- nia ciepła dostarczanego do układu. Przykładowe wyni- ki obliczeń pokazano na rys. 3.
Przedstawione na rys. 3 wyniki symulacji pokazują, że wzrost wielkości qd prowadzi do nasilanie się oscylacji wielkości opisujących rozpatrywany układ. Na rys. 3a pokazano wyniki symulacji przy qd = 85 kWm-2.Ampli- tuda oscylacji w kolejnych iteracjach zmniejsza się co prowadzi do zaniku oscylacji – w układzie ustala się stan ustalony. Na rys. 3b pokazano wyniki symulacji dla qd = 95 kWm-2.Amplituda oscylacji parametrów opisu- jących rozpatrywany układ zwiększa się w kolejnych iteracjach. Prowadzi to do pojawienia się w układzie niegasnących oscylacji parametrów opisujących wymianę ciepła i masy.
Rys. 3. Wpływ wartości qd na charakter oscylacji: G, qu. a) qd = 85 kWm-2 b) qd = 95 kWm-2.
Mechanizm pojawienia się w modelu oscylacji parame- trów opisujących układ jest następujący. Występujące w mikrokanale wrzenie powoduje zatykanie się mikrokana- łu powstającymi pęcherzami parowymi, prowadzi to do wzrostu spadku ciśnienia na długości mikrokanału.
Odblokowanie się mikrokanału następuje, gdy ciśnienie cieczy na wlocie do mikrokanału jest wystarczająco wysokie do zwiększenia przepływu cieczy do wartości wystarczającej do zmniejszenia intensywności wrzenia.
Jednak zmniejszenie intensywności wrzenia w mikroka- nale wywołuje wzrost temperatury powierzchni grzejnej, który prowadzi do wzrostu intensywności wrzenia - cykl oscylacji powtarza się. Częstotliwość oscylacji zależy od warunków przepływu ciepła i masy.
Gdy gęstość strumień dostarczanego do układu ciepła jest mała, wówczas zmiany warunków przepływu towa- rzyszące zmianom temperatury powierzchni grzejnej są wystarczające do wygaszenia oscylacji temperatury powierzchni grzejnej. W kolejnych iteracjach amplituda oscylacji maleje i w układzie ustala się stacjonarny stan wymiany ciepła i masy (rys. 3a). Jednak gdy strumień dostarczanego do układu ciepła zwiększa się, wówczas amplituda oscylacji temperatury również się zwiększa.
Takie oscylacje temperatury powierzchni grzejnej wy- magają większych zmian przepływu w kanale. Zmiany
przepływu w modelu opisane są równaniem (9). Wystę- pujące w równaniu (9) współczynniki oraz założenie o liniowości zmian ciśnienia na długości mikrokanału powoduje, że uzyskana z równania (9) reakcja układu na oscylacje temperatury powierzchni grzejnej jest niewystarczająca i układ traci stabilność (rys. 3b).
W zaproponowanym modelu prowadzi do pojawienia się niegasnących oscylacji.
3. PODSUMOWANIE
W pracy analizowano wymianę ciepła i masy w układzie zawierającym mikrokanał. Przedstawiono model pozwa- lający symulować oscylacje: natężenia przepływu, ci- śnienia, temperatury powierzchni grzejnej oraz natężenia strumienia ciepła pobieranego przez wrzącą ciecz. Do opisu właściwości przepływu mieszaniny dwufazowej zastosowano jednowymiarowy model przepływu typu –
‘homogeneous flow model’. Model przepływu ciepła w ściance mikrokanału oparty został na jednowymiaro- wym równaniu przewodnictwa ciepła z uwzględnieniem zmiennego w czasie strumienia ciepła pobieranego przez wrzącą ciecz.
Otrzymane wyniki symulacji pokazują, że pojawienie się oscylacji parametrów opisujących wymianę ciepła i masy w mikrokanale określona wzajemna interakcja
pomiędzy układem zasilającym wymiennik ciepła w ciecz a układem dostarczającym ciepło do mikrokana- łu. O stabilności takiego układu decyduje jego charakte- rystyka. W praktyce oznacza to, iż poprzez odpowiednie zmiany w konstrukcji układu zasilającego wymiennik w ciecz można doprowadzić do usunięcia oscylacji parametrów opisujących wymianę ciepła i masy.
Zaproponowany model pozwala określić graniczne warunki stabilnej pracy układu. Jednak uzyskane w ten sposób wartości odpowiadają układowi w którym wy- miana ciepła opisana jest korelacjami zastosowanymi
w prezentowanej pracy. Zastosowanie innych korelacji opracowanych przez innych badaczy zmieni wartości krytycznych parametrów stabilności. Dlatego uzyskany wynik można jedynie traktować jako jakościowy – opisujący oddziaływanie pomiędzy poszczególnymi mechanizmami wymiany ciepła i masy w układach z mikrokanałami.
Uzyskany w ramach zaproponowanego w pracy modelu scenariusz pojawiania się oscylacji parametrów opisują- cych wymianę ciepła i masy obserwowany był w bada- niach eksperymentalnych opisanych w pracy [4, 5, 13].
Literatura
1. Ledinegg M.: Instability of flow during natural and forced circulation. “Die Waerme” 1938, 61, p. 891–898.
2. Stenning A.H.: Instabilities in the flow of a boiling liquid. “Journal of Basic Engineering” 1964, 86, p. 213–228.
3. Kakac S.,Bon B.: A review of two-phase flow dynamic instabilities in tube boiling systems. “International Journal of Heat and Mass Transfer” 2008, Vol.: 51 Issue: 3-4.
4. Mosdorf R., Ping Cheng, Wu H.Y., Shoji M.: Non-linear analyses of flow boiling in microchannels. “International Journal of Heat and Mass Transfer” 2005, Vol. 48, p. 4667-4683.
5. Brutin D., Tadrist L.: Pressure drop and heat transfer analysis of flow boiling in a minichannel: influence of the inlet condition on two-phase flow stability. “International Journal of Heat and Mass Transfer” 2004, Vol. 47, p.
2365–2377.
6. Thome J. R.: Engineering Data Book III , Wolverine Tube, Inc ,2004.
7. Wu H.Y., Cheng P.: Visualization and measurements of periodic boiling in silicon microchannels. “International Journal of Heat Mass Transfer” 2003, Vol. 46, p. 2603–2614.
8. Wu H.Y., Cheng P.: Liquid/two-phase/vapor alternating flow during boiling in microchannels at high heat flux.
“Int. Commun. Heat Mass Transfer” 2003, Vol. 30, p. 295–302.
9. Wu H.Y., Cheng P.: Boiling instability in parallel silicon microchannels at different heat flux. “International Journal of Heat and Mass Transfer” 2004, Vol. 47, p. 3631–3641.
10. Warrier G.R., Dhir V.K.: Momoda L. A.: Heat transfer and pressure drop in narrow rectangular channels. “In- ternational Journal of Heat and Mass Transfer” 2008, Vol. 51, p. 399–433.
11. Cicchitti A., Lombaradi C,. Silversti M., Soldaini G., Zavattarlli R.: Two-phase cooling experiments – pressure drop heat transfer burnout measurements. “Energia Nucleare” 1960, 7, p. 407–425.
12. Carey V.P.: Liquid-vapor phase-change phenomena: an introduction to the thermophysics of vaporization and condensation processes in heat transfer equipment. Hemisphere Publishing Corp., Washington 1992.
13. Wang G., Cheng P., Wu H.: Unstable and stable flow boiling in parallel microchannels and in a single microchannel. “International Journal of Heat and Mass Transfer” 2007, Vol 50, p. 4297-4310.
