YATAY BORU ÜZERİNDE LAMİNER FİLM YOĞUŞMASININ NÜMERİK ANALİZİ
K. ALDAŞ ve K. ALTINIŞIK
Niğde Üniversitesi. Aksaray Mühendislik Fakültesi, Makine Mühendisliği Bölümü, 68100, Aksaray
ÖZET
Bu makalede; akışı yerçekimi ile sağlanan film yoğuşmasının, bir sayısal analizi yapılmıştır. Analizde, hidrodinamik ve termal sınır tabaka denklemleri kullanılmıştır. Momentum denkleminde bulunan atalet terimleri ile enerji denkleminde bulunan konveksiyon terimlerinin yoğuşmaya etkisi incelenmiştir. Su buharının 20 mm çapında bir yatay boru üzerinde yoğuşmasından oluşan sıvı film kalınlıklarının tüp yüzeyindeki dağılımı hesaplanmış ve Nusselt hipotezinden hesaplanan sıvı film kalınlığı ile karşılaştırılması yapılmıştır.
Anahtar kelimeler: Yatay boru, film, yoğuşma, laminer akış
NUMERICAL ANALYSIS OF LAMINAR FILM CONDENSATION ON A HORIZONTAL TUBE
ABSTRACT
In this paper a numerical analysis of gravity induced laminer film condensation is conducted. In the analysis hydrodynamic and thermal boundary layer equations are used. Effects of inertia terms in momentum equation and convective terms in energy equation are examined. The distribution of the liquid thickness on the tube surface caused by the condensation of steam on a 20 mm diameter tube is determined and compared with obtained from Nusselt hypothesis.
.
Keywords: Horizontal tube, film, condensation, laminar flow
1. GİRİŞ
Film yoğuşmasının ilk analizi Nusselt [1] tarafından yapılmıştır. Bu analizde, sıvı film içindeki sıcaklık dağılımı lineer kabul edilmiş ve yer çekim kuvveti ile viskoz kuvvetleri dengelenerek sıvı tabakanın kalınlığını veren bir denklem elde edilmiştir. Sparrow ve Gregg[2], Nusselt’in analizinde gözardı edilen atalet kuvvetlerini de hesaba katan bir analiz yapmışlardır. Aynı analizde, aşırı kızdırılmış buhardan sıvı film tabakasına olan konveksiyon enerjisi de hesaba katılmıştır. Sparrow ve Gregg’in analizleri, sıvı film içindeki sıcaklık dağılımının hemen hemen lineer olduğunu göstermektedir. Sıvı Pr sayısı 1 ve 100 arasında ise ve (Ts-Tw)<10 olduğu sürece, atalet terimlerinin önemsiz olduğunu ifade etmektedirler. Koh ve arkadaşları[3]’nın analizinde sıvı tabaka ile çevresindeki gaz ortam arasında yüzey sürtünmesinden dolayı meydana gelen momentum alışverişinin, sıvı tabakanın hareketine olan etkisi incelenmiştir. Normal endüstriyel uygulamalarda (soğutma tekniği ve kondenser tasarımı gibi) sıvı tabaka ile çevresi arasındaki gaz ortamı arasında meydana gelen momentum alışverişinin önemsiz olduğunu belirtmişlerdir.
Sparrow ve Eckert[4], kızgın buharın ve buhar içindeki yoğuşmaz gazların ısı taşınım katsayısı üzerindeki etkisini incelemişlerdir. Kızgınlık derecesinin sebep olduğu ısı transferinin su ve Freon tipi soğutucular için çok
önemsiz olmakla birlikte, yoğuşmaz gazların ısı transferini küçülttüğü görülmektedir. Sparrow ve Lin[5]’in yoğuşmaz gazların etkisi üzerine yaptıkları analizler, yoğuşmaz gazların ısı transferini büyük ölçüde düşürdüğünü göstermektedir. Shekriladze ve Gomelauri[6], buhar süpürme etkisi altındaki yoğuşmayı ilk olarak analiz etmişlerdir. Yapılan analizde, yoğuşmanın sebep olduğu yoğuşma yüzeyine dik hız bileşeninin buhar tabakasındaki hız dağılımına önemli ölçüde etki ettiği, ara yüzeydeki hız gradyantını büyüttüğü ve neticede buhar süpürmesinin arttığı gösterilmiştir. Fujii ve Uehara[7], Shekriladze tarafından ele alınan buhar süpürmesi olayını daha gerçekçi şartlarda incelemişlerdir. Gaddis[8], tarafından yapılan analizde hem sıvı tabaka hareketi hemde buhar sıvı tabakasının hareketi, komple sınır tabaka denklemleriyle ele alınmıştır. Bu analiz, o zamana kadar yapılan en kapsamlı incelemedir. Karabulut ve Ataer[9], hem yoğuşmuş tabakanın hareketini hem de buhar sıvı tabakanın hareketlerini komple sınır tabaka denklemleri ile simüle etmişlerdir. Önceki araştırmalardan farklı olarak denklemlerin çözümlenmesinde, sonlu farklar metodu kullanmışlardır. Muhtelif şartlarda ısı taşınım katsayısının dağılımını vermişlerdir. Srzic ve arkadaşları[10], eğik konumlu izotermal bir pleyit üzerindeki buhar süpürmeli yoğuşma olayını, yoğuşmaz gazların varlığını dikkate alarak incelemişlerdir.
2. FİZİKİ MODEL VE MATEMATİK MODEL
Yerçekimi kaynaklı olayın fiziki mekanizması Şekil 1 de görülmektedir. İçerisinden soğutucu su geçen borunun çevresinde bulunan doymuş durgun buhar, borunun cidarlarında yoğuşmaktadır. Yoğuşan madde belirli bir kalınlığa ulaştıktan sonra, laminer rejimde yerçekimi etkisiyle aşağıya doğru akmaktadır. Borunun cidarlarında bir film tabakası oluşturan yoğuşmuş madde, borunun cidar kalınlığı ile birlikte ısıl direnç oluşturmaktadır. Yoğuşma ısı transferinde araştırılan husus, sıvı film ısı direncidir. Isı transferinin analizi için sıvı tabaka içerisindeki sıcaklık dağılımının ve sıvı tabaka kalınlığının belirlenmesi gerekir. Sıvı tabaka hareket halinde olduğu için, kalınlığının belirlenmesi yoğuşan miktarı ile akan miktarın dengelenmesi gerekir.
i c
o m m
m& = & +& (1) Eşitlik 1 de birim zamanda
m &
i kütle balansı için seçilen 1 m derinliğindeki elemana girer kütleyi, m&celemanın serbest yüzeyinde yoğuşan madde miktarını, m&o ise elemandan çıkan kütleyi göstermektedir. m&o yerine,∆x
dx m m d m
i i
o
+
= &
&
& (2)
yazılabilir. Buna göre 1 nolu denklemdeki kütle balansı;
mc
dx x m d& &
=
∆ (3)
şeklini alır. 3 nolu eşitlikte bulunan m&c, (1.∆x) genişliğindeki yüzeyde birim zamanda meydana gelen yoğuşma miktarı olup;
h ∆x
∆x q
"
m m
fg c
&
& = = (4)
şeklinde ifade edilir. 2 nolu eşitlikte bulunan m&yoğuşmuş maddenin akış hızı cinsinden
Şekil 1 Fiziki model ve koordinat sistemi
∫
=
δ
0
udy
m& ρ (5)
şeklinde ifade edilir. 2, 3 ve 4 nolu eşitlikler birlikte değerlendirilirse, sıvı film kalınlığını belirleyen eşitlik, h ∆x
∆x q udy dx ρ
d
fg δ
0
&
=
∫
( 6)olarak belirlenir. Son eşitlikte bulunan q& yoğuşan maddenin sıvı film serbest yüzeyinde bıraktığı ısı olup, bunun bir kısmı akışkanın hareketinden dolayı, x yönünde taşınmaktadır. x yönünde taşınan ısı katı cidara geçen ısıya kıyasla çok küçük olduğundan ihmal edilir. Katı cidara geçen ısı Fourier yasasına göre
0 y
f dy
k dT q
=
=
& (7)
eşitliği ile hesaplanabilir. 7 nolu eşitlik 6 nolu eşitlikte yerine konursa
0 fg y
f y δ
o
y dy
dT h udy k dx ρ
d
=
=
=
=
∫
(8)olur.
δ
η= y (9)
Tw
mc
x X
Y
x+∆x g Ts
ro
v u δ mi
mo
Şekil 2 Dönüştürülmüş koordinat sistemi
şeklinde bir koordinat tanımlarsak integralin y=o sınırına karşılık η=0 , y=δ sınır şartına karşılık η=1 olur. Buna karşılık denklem 8 ise,
( )
0 fg η
f 1
0 dη
dT δ 1 h dη k uδ dx ρ
d
=
=
∫
(10)olur. İntegral sınırları x’ den bağımsız olduğu için d/dx operatörü integralin önüne geçebilir.
( )
η 0 fg
f 1
0 dη
dT δ 1 h dη k uδ dx ρ
d
=
=
∫
(11)Sonuç olarak,
η 0 1
0 fg
f 1
0 2
dη dT ρh udη k dx dδ dη δ x δ u
=
−
∂ +
∂
∫
∫
(12)elde edilir. Son eşitlikte bulunan d 0
dT
=
η η ın belirlenmesi için sıvı film içindeki sıcaklık dağılımının belirlenmesi gerekir. Bunun için termal sınır tabaka denklemi
2 2
y α T y v T x u T
∂
= ∂
∂ + ∂
∂
∂ (13)
kullanılabilir. Hem son eşitlik hem de 12 nolu eşitlik sıvı tabaka içindeki akış hızının belirlenmesini gerektirmektedir. Sıvı tabaka içindeki akış hızının belirlenmesi için hidrodinamik sınır tabaka denklemleri
+
∂
= ∂
∂ + ∂
∂
∂
2 o 2
r sin x y g ν u y v u x
u u (14)
0
y v x
u =
∂ +∂
∂
∂ (15)
kullanılabilir. Eşitlik 12,13,14 ve 15 eş zamanlı çözülmesi gereken bir denklem takımı olup sınır şartları aşağıdaki gibi olacaktır.
x=0 , 0 dx
dδ= , 0
x T=
∂
∂ , u=0 (16)
y=0 , T=Tw , u=0 , v=0 (17)
y=δ , T=Ts , 0
y u=
∂
∂ (18)
Denklem 13,14 ve 15 in çözüm bölgesi sıvı filminin içerisidir. Fakat sıvı film kalınlığı x ile değiştiği için çözüm bölgesi düzgün değildir. Bu sebeple 13,14 ve 15 eşitliklerinde
δ
η= y koordinatı kullanılarak çözüm bölgesinin Şekil 2 deki gibi normalleştirilmesi gerekir. Yoğuşmayı yöneten denklemler
2
2
2 η
T δ
α η T δ v x u T
∂
= ∂
∂ + ∂
∂
∂ (19)
+
∂
= ∂
∂ + ∂
∂
∂
2 0 2
2 r
sin x g η
u δ
ν y u δ v x
u u (20)
0
η v δ 1 x
u =
∂ + ∂
∂
∂ (21)
şeklini alır. Eşitlik 13, 14, ve 15 in sınır şartlarını yeniden düzenlersek
x=0 , 0
dx
dδ= , 0,
x T=
∂
∂ u=0 (22)
η=0 , T=Tw, u=0, v=0 (23)
η=1 , T=Ts 0
η u=
∂
∂ (24)
olur. Denklemlerin sonlu farklar şekli ekte verilmiştir. Nusselt sayısı ise aşağıdaki denklem ile hesaplanmıştır.
( ) (
∞)
=
∞
=
−
− =
= T T
δ D dη dT
T T dy D dT Nu
s η 0 s
0
y (25)
3. BULGULAR VE TARTIŞMA
Şekil 3’de ∆T=10 °C için belirlenen hız dağılımının, muhtelif açılardaki profili görülmektedir. x=0 civarında serbest yüzeydeki akış hızı yeterince küçükken x’in ileri değerlerinde, serbest yüzeydeki akışın hızlandığı görülmektedir, Burada; yoğuşma olayının analizinde momentum denklemlerinde bulunan atalet terimlerinin özellikle özgül ağırlığı fazla olan akışkanlarda önem kazandığı görülmektedir, Momentum denklemlerinde bulunan atalet terimleri ihmal edilirse, borunun alt tarafında bulunan sıvı tabakanın kalınlığı olması gerektiğinden daha hızlı büyüyecektir. Gerçekte ise, borunun üst yüzeyinde kazandığı hızla, borudan daha yüksek bir hızla ayrılacak ve δ kalınlığı aşırı büyümeyecektir.
Şekil 3 ∆T=10°C için değişik açılarda Eta-Hız profilleri
Şekil 4’de ∆T=10 °C için çeşitli açılarda belirlenen η istikametindeki sıcaklık profilleri birbirleriyle kıyaslanmaktadır. Profiller birbiri üzerine bindiği için değişik açı değerlerinde sıcaklık profilleri fazla değişmemektedir. Sıcaklık profillerinin lineer olması ve birbirine denk olması, akış istikametlerinde akışkanın hareketlerinden doğan enerji taşınmasının önemsiz olduğunu göstermektedir. Bu sebeple su buharının yoğuşma ısı transferi analizinde Nusselt hipotezinde yapılan lineer sıcaklık dağılımı kabulü kusursuz olmaktadır. ∆T=40
°C için yapılan incelemede aynı gidişatı göstermektedir.
Şekil 4 Değişik açılarda Eta boyunca sıcaklık dağılımı
Şekil 5’de Nusselt sayısının boru çevresindeki dağılımı muhtelif ∆T için verilmiştir. ∆T küçüldükçe Nusselt sayısı büyümektedir. Nusselt sayısının maksimum değeri borunun en üst noktasında görülmektedir. En üst noktasından uzaklaştıkça δ büyüdüğü için, ısıya karşı bir direnç oluşmaktadır ve Nu sayısı azalmaktadır.
0 0.2 0.4 0.6 0.8 1
25 30 35 40
S ıc aklık (T s -T w )
Eta
5 45 90 120 160
Şekil 5 Farklı ∆T’ler için Nusselt sayısı-Açı pozisyonu dağılımı
Borunun en alt noktasında Nu sayısı hesaplanamamakla birlikte sıfır olacağı görülmektedir. ∆T artarken q& ‘da arttığı için yoğuşan madde miktarı da artmakta, neticede borunun tepe noktasında δ, ∆T’ye bağlı olarak daha büyük değer almaktadır. Bu da Nu sayısının küçülmesine sebep olmaktadır.
Yoğuşma ısı transferi olayında, ∆T’yi büyük tutarak daha fazla ısı transferi gerçekleştirmek mümkün olmamaktadır. Bir optimizasyonla ∆T’yi arttırmanın nereye kadar faydalı olacağının incelenmesi gerekmektedir.
Şekil 6 Farklı ∆T’ler için Yoğuşma kalınlığı(δ)-Açı pozisyonu(φ) dağılımı
Şekil 6’da sıvı tabakanın kalınlığının boru çevresindeki dağılımı, muhtelif ∆T için verilmektedir. δ nın minimum değeri boru tepe noktasında görülmektedir. 90° civarında δ’nın küçülmesi beklenebilir. Çünkü 90°
civarında akışkanı hareket ettiren kuvvet maksimumdur. Fakat su buharı için yapılan analiz δ ‘nın sürekli arttığını göstermektedir. ∆T büyük olduğu zaman δ da bağıl olarak büyük olmaktadır ve q&’nun istenildiği gibi büyük
0 1 00 2 00 3 00 4 00 5 00 6 00
0 3 0 6 0 90 12 0 1 5 0 18 0
A çı p ozis yo nu (d ere c e)
Nusselt sayısı (Nu)
2 c 5 c 10 c 20 c 40 c
0 2 4 6 8 10 12 14 16 18 20
0 30 60 90 120 150 180
Açı pozisyonu (derece)
Yoğuşma kalınlığı (*E-5) (mm)
2 c 5 c 10 c 20 c 40 c
Şekil 7 ∆T=10°C için muhtelif açı değerlerinde Nusselt hipotezi ile mevcut Nu sayısından elde edilen yoğuşma kalınlığının karşılaştırılması.
olmasını önlemektedir. Borunun alt yüzeyinde en alt noktaya doğru yaklaşan δ sonsuza doğru gitmekte ve 172°
nın daha ilerisindeki girid çizgilerinde, çözümleme işlemine devam edilebilmesi için 1/200 den daha küçük dη kullanmayı gerektirmektedir.
Şekil 7’de ∆T=10 °C için Nusselt hipotezi ve mevcut modelden çıkan δ kıyaslanmaktadır. Yegane fark 90°
den sonra Nusselt modelinde δ’nın daha hızlı büyümesidir. Bu atalet terimlerinin ihmalinden kaynaklanan bir farktır.
4. SONUÇ
1. Su buharının yoğuşması olayında sıvı tabaka içerisindeki sıcaklık dağılımı lineer ve enerjinin konveksiyonla akışkan tarafından taşınımı önemsizdir.
2. Atalet terimlerinin 90° den sonra az miktarda fark oluşturduğu görülmektedir.
3. Yoğuşma olayında ∆T yı yeterince büyük tutarak ısı akısını yükseltmek mümkün değildir.
4. Elde edilen sonuçlar Nusselt hipotezinden elde edilen sonuçlara çok yakındır.
SEMBOLLER Cp Sabit basınçta özgül ısı (J/kgK)
D Boru çapı (m)
g Yerçekim ivmesi (m/s2) hfg Buharlaşma entalpisi (J/kg) k Isı iletim katsayısı(W/mK) Nu Nusselt sayısı D(dT/dy)w/(Ts-Tw) m& Kütle debisi (kg/s)
T Sıcaklık (°C) Pr Prandtl sayısı (µCp/k) q& Isı akısı (W/m2)
ro Boru yarı çapı (m) u x doğrultudaki hız (m/s) v y doğrultudaki hız (m/s) α Isıl yayılım (m2/s) δ Film kalınlığı (m)
η Boyutsuz koordinat(y/δ)(Eta) µ Dinamik viskozite (kg/m s) ν Kinematik viskozite (m2/s) ρ Yoğunluk (kg/m3)
0 5 10 15 20 25 30
0 30 60 90 120 150 180
Açı pozisyonu(derece)
Yoğuşma kalınlığı(*E-5)(mm)
Mevcut Nu sayısı Nusselt hip.
KAYNAKLAR
1. Nusselt, W. ( 1916.) ”Die Oberflachen-Kondensation des Wasserdampfes”,VDI Z.,vol.60, pp. 541-546 and 569-575.
2. Sparrow, E.M.,and Gregg, J.L. ( 1959) ”A Boundry-Layer Treatment of Laminar-Film Condensation“, J.of Heat Transfer, 81C,13-18.
3. Koh, C.Y.,Sparrow,E.M.,and Harnett, J.P.(1961) ”The Two Phase Boundary Layer in Laminar Film Condensation”,Int. J.Heat Mass Transfer 2, 69-82.
4. Sparrow, E.M.,and Eckert (1961) ”Effects of Superheated Vapor and Noncondensable Gases on Laminar Film Condensation”,AIChE Journal, vol.7,pp.473-477.
5. Sparrow, E.M., and Lin, S.H.(1964) ”Condensation Heat Transfer in the Presence of a Noncondensable gas”,Transaction of the ASME,pp.430-436.
6. Shekriladze, I.G. and Gomelauri, V.I.(1966) ”Theortical Stady of Laminer Film Condensation of Flowing Vapour”, Int. J. Heat Mass Transfer vol.9, 581-591.
7. Fujii, T. and Uehara, H.(1972) “Laminar Filwise Condensation on a Vertical Surface”, Int. J. Heat Mass Transfer.,vol.15,pp.217-233.
8. Gaddis, E. S.(1979) “Solution of the Two Phase Boundary-Layer Equations for Laminar Film Condensation of Vapour Flowing Perpendicular to a Horizontal Cylinder”,Internationl Jornal of Heat Mass Transfer ,vol.22,pp.371-382.
9. Karabulut, H. ve Ataer, Ö.E.(1996) ”Numerical Analysis of Laminar Film-Wise Condensation”,Int J.
Refri. Vol.19, No. 2, pp. 117-123.
10. Srzic,V.,Soliman, M. and Ormiston, S.J.(1999) ”Analysis of laminar mxed-convection condensation on isothermal plates using the full boundary –layer equation mixtures of a vapor and a lighter gas”, Int, J.
Heat and Mass Transfer,vol.42 pp.685-695.
11. Curtıs F, Gerald,”Appliied Numerical Analysis”, Addison-Wesley Publishing Company , Londen
EKLER
Kullanılmakta olan non-lineer denklemler için y istikametinde zorunlu olarak bir iteratif çözümleme işlemi uygulanmıştır. Sonlu fark eşitlikleri, sayısal yöntem, ve sayısal çözümleme işlemlerinde Over Relaxation[11]
metodu kullanılmıştır. Eşitliklerin sonlu fark eşitlikleri aşağıda verilmiştir.
Momentum Denklemi:
+
∂
= ∂
∂ + ∂
∂
∂
2 0 2
2 r
sin x u g u
v x u u
η δ
ν η δ
[ ]
0r sin x g U U
∆η δ U ν
∆ηδ U v
∆x u
∆η δ
2ν δ
∆η v
∆x U u RESU
0 1
j i, 1 j 2 i,
1 j i, j
1, 2 i
j 2
i, − + − + − =
− +
= − + + −
− +
−
=
2 2 j
i, j i,
∆η δ
2ν δ
∆η v
∆x u U RESU
U (26)
Süreklilik Denklemi:
η 0 v δ 1 x
u =
∂ + ∂
∂
∂
[ ] [
V V]
0η δ
∆ U 1
∆x U
RESV= 1 i,j− i,j−1 + i,j− i,j−1 =
[
RESV]
V δ∆η
Vi,j= i,j− (27)
Enerji Denklemi:
T 0 η α T δ v x
u T 2
2
∂ =
− ∂
∂ + ∂
∂
∂
η
[ ] [ ] [
T 2T T]
0∆η δ 1 α j Ti, j δ∆η Ti, T v
∆x T
REST= u i,j− i−1,j + − − − 2 2 i,j−1− i,j+ i,j−1 =
+ +
−
=
2 2 j
i, j i,
∆η δ
2α δ∆η
v
∆x u T REST
T (28)
( )
2 1
1
0
1
0 1
0 f fg 0
1 f 2n 2
xd x u 2 ud
d dT h x k 2 ud
∂
∆ ∂ +
∆ +
=
∫ ∫
∫
=
−
η η
η η ρ
δ
δ η (29)
olmaktadır.
Çözümleme işlemi x=0 da bulunan birinci grid çizgisinden başlayıp x istikametinde çizgi çizgi devam ettirilmiştir. İlk çizgi u=0 olarak belirlenmektedir. Bu çizgi üzerinde T ve δ bilinmemektedir. Bu sebeple δ ve T’nın bir tahmini yapılarak ikinci grid çizgisi üzerindeki u, v,T ve δ belirlenebilir. Son belirlenen u,ve T dağılımı 29 nolu eşitlikte kullanılarak δ hesaplanır. Bu hesaplama δ’nın hesabına başlarken kullandığımız atma değerden daha doğrudur. Belirtilen δ kullanılarak u, ve T ikinci grid çizgisi üzerinde yeniden belirlenir. Verilen sınır şartının gereği olarak 1. ve 2. Girid çizgisi üzerindeki T dağılımı birbirinin aynısı olacaktır. Bu şekilde δ hesabı birkaç kez tekrarlandıktan sonra, aynı işlemler üçüncü ve müteakip girid çizgilerine uygulanarak hesaba devam edilir.
