矩形多孔質層の非定常凍結熱伝達
佐々木 章
TransientFreezingHeatTransferof aWater‑Saturated PorousMedia inaVerticalRectangularCavity
AkiraSAsAKI
(昭和63年10月31日受理)
Thepurposeofthepresentstudyistoexamineanalyticallyandexperimentallythe effectsofthebeaddiameterandthephysicalpropertyofporousmediaontransientfreezing heat transferofawatersaturatedporousmedia・Thephysicalmodelconsideredinthis studyisatwo‑dimensionalverticalcavity・Thetemperaturesofbothhotandcoldwalls whichisverticallyoppositeareuniform,respectively・ Theequationofmomentumtakes acmunt㎡Fordlheimm'sextensionastheresistancetonowintheporousmedia・Thegoverning equationsaretransformedbytheboundaryfixingmethod. TheSORmethodisutilizedto solvenumericallytheequations・Theexperimentsareperformedinarectangulartestcell thathaveinsidedimensionsoflOcminheightanddepthand5cminwidth,Threedifferentsize glass,ironandcopperbeadswereusedastheporousmanialinthisstudy.Thetemperature ofthecoldwallwas‑10℃andoneofthehotwallwas12℃.Comparisonsoftheanalytical resultswiththeexperimentalonesaremade,andtheeffectsofthebeaddiameterandthe physicalpropertyofporousmediahavebeendiscussedfortransientfreezingheattransfer.
1. 緒
言巨
時間帯に冷房に利用する氷蓄熱空調システムが利用 されている。氷蓄熱は状態変化を伴う潜熱蓄熱であ るため,従来の温度変化を利用した顕熱蓄熱である 水蓄熱に対し蓄熱量の増大だけでなく,省スペース をはかることができる。氷蓄熱槽を開発する上で,
考えるべき事は,熱負荷に見合った蓄熱を行い,日 中の追加運転を最小限に抑え,料金の安い夜間電力 のメリットを最大限生かす制御ができることである。
すなわち,製氷,解氷時の熱伝達特性の制御である。
しかし, この蓄熱槽の性能は利用する蓄熱材料すな わち水,氷の熱的性質によって限定されてしまい,
その利用目的に応じた伝熱特性を得ることが困難で ある。 しかし,水と固体粒子よりなる多孔質層を利 用した蓄冷槽の場合,その利用目的に応じて固体粒 子の種類,粒子量を変えることにより幅広くその定 常蓄熱特性を変ええる可能性がある(2)。したがって,
凍結を伴う場合の粒子層の伝熱特性を知ることは,
自然凍結に伴う凍害防止の基礎となるだけでなく,
蓄冷槽の開発設計の基礎資料となる。多孔質層の凍 結に関する研究(3) 伯)はこれまでにもなされている 含水した多孔質層の凍結現象は,寒冷地での土壌
の自然凍結に代表される他,工業的にも軟弱な地盤 の固化を目的とした地盤凍結工法,食品の凍結保存
・加工などに見受けられる。土壌の自然凍結は,地 下埋設物すなわち水道管や下水道管の凍結閉塞を引 き起こす原因となるだけでなく,地下建造物の断熱 の問題とも密接に関連する。これらの凍結現象によ る.弊害を避けることができれば,社会生活上大きな 貢献となる。一方,近年エネルギーの有効利用とい う観点から寒冷地においても地熱熱源ヒートポンプ の利用が活発に試みられている。この種のヒートポ ンプの性能は,地下埋設熱交換器まわりの土壌の凍 結や融解に影響されることが報告されている(')。ま た,同様の立場から土壌を潜熱蓄熱槽として利用し,
冬期の冷熱を夏期に有効に活用しようとする試みが なされている。また,都市部においては,昼と夜で 利用する電力に相当量の差がみられることから,深 夜電力により氷蓄熱を行い, 日中の熱負荷の大きい
−13−
矩形多孔質層の非定常凍結熱伝達
が,実験範囲が限られているうえ,凍結熱伝達特性 に関してはほとんど検討がなされておらず,十分な 知見が得られているとは言えない。
以上の観点から,本研究は,最も基本的なモデル である矩形多孔質層の凍結熱伝達現象に及ぼす粒子 径,および粒子物性の影響を数値解析ならびに実験 により検討しようとするものである。
2. 数値解析
図1に座標系を示す。高さH,幅Wの矩形空間を 考え,上下面を断熱壁,左右の垂直面をそれぞれ加 熱壁,冷却壁とする。従来,多孔質層内の流体の流 れには,ダルシーの法則を適用することが多い。そ の結果,高レーレー数領域において計算結果と実験 結果とに差異が生じると指摘されている。また,ダ ルシーの法則では、壁面ですべりを許すことになる。
本解析では,壁面でのすべりなしの条件を得るた め粘性の影響(Brinkmanモデル(7))を, また高レ
ーレー数領域への適用が可能となるよう速度の二乗 に比例した慣性の影響(Forchheimerモデル(8))を
取り入れた運動方程式を採用し謝解析に当たり次
の仮定を導入する。 (1)流体の流れは, 2次元・非圧 縮・層流である。 (2)熱物性値は,浮力の項に含まれ る密度を除き一定である。 (3)壁面近傍の空隙率の不 均一は考えない。多孔質層は均一である。
以上の仮定の基で,基礎方程式は次のようになる。
液相(未凍結粒子層,以後液相と呼ぶ)領域での 連続の式,運動方程式,およびエネルギ方程式は,
記 号
慣性係数 比熱
ダルシー数=k/W2 粒子径
液相厚さ(凍結界面位置)
重力加速度 高さ 熱伝達率 浸透率 潜熱
ヌセルト数=hf/ス 圧力
温度 凍結温度 時間
ダルシアン速度 幅
座標 空隙率
無次元温度=(TI‑To)/(Th‑To) 熱伝導率
粘性係数 密度
無次元時間(フーリニ数) =(・端,
流れ関数 加熱壁 冷却壁 凍結界面
:::..::::::::::::::::..:..:..::・・字:::::::●●WWFVV︾
auO99.︒・
CCDdfgH・hkLNPTTtuWXE6i〃βて〃HCI添CeflhmS抑一触
一一j伽諏
V+
伽 一 u 献
/11画︑
1uljノpl夢0ノーI︑
一一十
伽一帥伽一肌
十︑1jノ伽一触〃催い
(1)
1 12
偲仰両十一u一一
糾紛
十十ん伴い伽一伽I卜1V巾八
牌紳 哨㈲
0
Th Tb
冷却壁 等価 凍結界面
液相(未凍結層)または局所 加熱壁
平均
固相(凍結層)
倖聖一ョ
図1 座標系
e儂寺劉佳・斜 3
+
(。,).!評十(cp)1(。語報課)
(謬半箒)
=スel
となる。固相(凍結粒子層,以後固相と呼ぶ)領域 でのエネルギ方程式は,
(cp).。f 。(帯寺辮) ⑮)
となる。凍結界面でのエネルギバランス式は,
)( )
(伽。。詩‑卿。,鶚叶儲‑・,。Lg= {6)
となる。初期条件,境界条件は次のようになる。
t=0 ;u=v=0 ,TI=Ts=Th
l ㈹
y=0 ;u=v=0 ,T,=Ts
y=f ;u=v=0 ,TI=Ts=0
y=W;Ts=Tc
x‑0,H;U、‑0, ;¥‑g¥‑'
試験部①は高さH100m,奥行き100mg,幅W50mで あり,厚さ10mのアクリル樹脂板で作成した。試験 部正面には,温度分布の可視化および凍結界面観察 のため厚さ1mのアクリル板を用いた。加熱および 冷却壁⑤⑥はいづれも100×100m,厚さ5mmの銅板 を用いた。加熱器として,短冊状に加工した厚さ30
〃mのステンレス箔を用いた。主加熱器⑨背面には 厚さ10mのベークライト板⑦を介して補償用加熱器
⑩を設けた。等温壁を得るために,主加熱器および 補償用加熱器は高さ方向に5分割し,各々独立に制 御できるようにした。冷却部は,冷却面⑥ 4分割 された冷却室②および多数の穴を開けたブライン噴 き付け用銅管③より構成されている。
加熱,冷却壁温度および主・補償用加熱器の温度 は直径0.3mのアルメルークロメル熱電対で測定し た。粒子層内の温度分布は,赤外線カメラを用いて 可視化した。充填粒子として,平均直径2mm(空隙 率e=0.36), 5m(e=0.42),および16nm(E
=0.48)の3種類のガラス球,平均直径2mの鉄球 (E=0.39)および銅球(E=0.39)を用いた。
実験条件は,冷却壁温度'It=‑10℃,加熱壁温度Tin
=12℃である。実験装置は周囲温度の影響をできる だけ軽減するために,実験中,厚さ50mのスタイロ フォーム断熱材でおおった。
加熱壁面,凍結界面での熱伝達率は,代表温度と して加熱壁温度および凍結界面温度を用いて次式の ように定義した。
(TM==。)Hf.¥スel
hmh ==‑
y=0d× (8)
征m4:+。)H〃い(簿)
hmf =‑
割y壽fd× (9)
式(2), (3)における浸透率k,慣性係数Cは,多孔 質層の構造によって変化するものであり,本解析で はErgun(9)により求められた結果を採用した。また,
各相の等価熱伝道率はKunii‑Smithの式('の上り算 定した。上記基礎方程式は流れ関数,渦度を導入し た後,境界固定法('')を用いて座標変換を行った。
基礎方程式の差分化はコントロールボリュウム法('2)に より行い,差分式はSOR法を用いて解いた。なお,
対流項には風上差分を適用した。,
〃ヨヨヨヨb
。⑩⑦⑨⑤⑥
①②③④⑤
⑥Copper plqte
⑦BqkeliteplOte
⑥Aluminumplqte
⑨Moinheoter
⑩Guqrdheqter TestsectiOn
Brineroom Brinepipe(in) Brinepipe(out) Copperplqte Thermocouples 3. 実験装置および実験方法
実験装置の概略を図2に示す。本実験装置は、加
熱部,試験部,および冷却部より構成されている。 図2実験装置
̲L正
一口デー]
│ │cgl‐‐q■
①
pーーq
■■一一 ■口一○一一一︒︽一一一一︒■重ロ■
■一一
■一一■
ll l
(<−可
rI−III
、
I
,
1 1、
②
1
0
○
○
ノわ
−15−
矩形多孔質層の非定常凍結熱伝達
I
=
e =1.0
ー
I
言.エ 三・○ 芸・工 三・︒
I
(b) (c)
(a)
図3 等温線および流れ関数(1) (Th=12℃,Tt=‑10℃,d=16m,t=5min)
e=1.0
A6=02
﹃三・工 三・○
三・工 三○
LL
ー
e=
0.0
(b) (c)
(a)
図3 等温線および流れ関数(2) (TIn=12℃,Tb=‑10℃,d=16m,t=lOmin)
匹・一俳虻
I
LL:
ー
蕊
e=1.0 Ae=0.2三.エ 三・︒
﹃三・○三・工
L昌
一
L昌
一
6 6=00
(b) (c) (a)
図3 等温線および流れ関数(3) (T1,=12℃、Tb=‑10℃,d=16m,t=30min)
凍結界面および加熱壁での局所熱伝達率h! f/hmf, h1h/hmhの計算結果を図4に示す。図はTI,=12℃,
Tc=‑10℃, t=30minでの結果でガラス粒子の場合 にたいするものである。d=2mの場合, hlf/hmf,h!h /hmh, いずれも高さ方向への変化はほとんど認めら れず,熱移動は主に伝導伝熱で行われていると言え る。h,h/hmhは粒子径が増大すると,右下がりの勾配 を持ち, d=5nunとd=16mでは似た分布となっている。
h,f/hmfはd=5mでは上部から下部Iご向かって単調に 減少するのに対し, d=16mではx/H=0.3付近で最 小値を示す。 これは,図3(c)から判断できるように,
下部に反時計方向まわりの流れが生じるためと考え られる。
図5は,平均熱伝達率hnmにおよぼす粒子径の影響 を示したものである。図より粒子径の増大とともに,
hmhの時間に対する増加割合が増大していることがわ かる。また, d==5nnn, 2nnnの場合, hmhは時間に対し増 加し続ける。 これは,液相領域の熱移動が主に伝導 伝熱により支配されるためと考えられる。すなわち 時間の経過にともない,固相領域に較べ等価熱伝導 率の小さい液相領域が減少するためと思われる。そ れにたいし, d=16mの場合はt=0.4hour付近で一 定値に達している。図6は熱物性の異なる粒子を用 4. 結果および考察
加熱壁温Th=12℃,冷却壁温Tt=‑10。Cで,粒子径 d=16mの場合の経時変化を図3(1)〜(3)に示す。 (a)は 温度分布の可視化写真で黒く写っている領域が0℃
以下の温度領域(凍結層)を表している。 (b)、c)は 計算により得られた等温線,流線である。いづれの 場合も,温度分布および凍結界面形状の計算結果と実 験値は良く一致していることがわかる。 t=5minで は,流れ関数(c)より加熱面で上昇,凍結界面で下降 する大きな循環流れが形成されているのがわかる。
さらに,下部の凍結界面近傍に,密度逆転による反 時計方向まわりの渦が発生しているのが認められる。
温度分布(a), (b)には,時計方向まわりの流れに伴い 熱が拡散している様子が示されてい、る。時間の経過 につれ,時計方向まわり大きな循環流のため,凍結 界面上部での熱伝達率が増大し, この領域の凍結の 進行が妨げられている。これは,等温線の間隔が上 部で粗く ,下部で密になっていくことからも判断で きる。また,反時計方向まわりの流れの領域が広が るとともに,流れ関数も増大している。方向の異な る流れが分離する付近で等温線の間隔が広くなって いる。
−17−
矩形多孔質層の非定常凍結熱伝達
いた場合のhmhとtの関係を示す。粒子の物性の影響 が顕著に現れていることがわかる。すなわち,粒子 の熱伝導率が大きいほどhmhの時間に対する増加割合 が大きくなり,短時間で一定値に到達している。計 算結果は,ガラスにおいては実験結果と良く一致す るが,粒子の熱伝導率が大きくなると実験結果より 大きな値を示す。特に銅粒子の場合には差異が大き い。これは,銅粒子と水の組合せの場合は両者の熱 物性値の差が極めて大きいためと考えられる。すな わち,液相,=固相の等価熱伝導率の推定の難しさ,
ならびに本解析で適用した多孔質層を熱的に均質な 物質とした仮定が成立しなくなるためではないかと 考えられる。
凍結量を潜熱蓄熱量Qs/Qで表し,その経時変化を ガラス粒子径をパラメータとして計算した結果を示 したのが図7である。Qsは凍結層の潜熱量,Qは試 料全体が凍結したときの全潜熱量である。d=5,2nnnで はQs/Qにそれほど違いは見られない。それに対し,
d=16mの場合,対流による移動熱量が多くなること から,Qs/Qは減少している。また, hmの増加割合 が粒子径により異なるにもかかわらず,Qs/Qはt
<0.2hourで粒子径の影響をほとんど受けないこ とがわかる。これは,いずれの条件においても,冷 却開始直後の熱移動が伝導伝熱で支配されているた めと考えることができる。d=16mではt=0.2付近 から自然対流の影響が現れ始め,Qs/Qの増加割合 がd=2, 5mの場合に比べ,減少する。図8には,
粒子の物性がQs/Qに及ぼす効果を示してある。粒 子の熱伝導率の増大とともに,短時間で定常状態に 達するとともに,Qs/Qも増大していることがわか
る。一方,液相と固相の等価熱伝導率の比(スes/スeD が約2.0と等しいガラス,鉄粒子の場合には,QJQ
︒m︵Unm︑︒0mmm一■由■■U二t正此Ce℃R 舵吋帥珈 mtmm S w−
●一■■■■︽■■■■■唾附牛昨IXGEOの一一一
120
08ヱEヘ三舌と
40
0 2 3
t hour
図5平均熱伝達率の経時変化
(粒子径の影響)
嵐102
4
Th=12。C,Tb=‑10。C.d=2mm Experiment‑Prediction
oGIQss の Iron
●Copper /Copper
10
86ヱE一三 420二E二
津 Iron
① ① の
。̲。̲銘竺s≦
2 3
t hour
図6平均熱伝達率の経時変化
(粒子物性の影響)
4
1.0
Q8
GIossbeods,Th=12。C,Tb=‑10。C 1.0
08 06 04 0.2
帖鮎エー×
rPdlrllr
①く⑩①
02
0 1 2 3 4
hlf/hmf ,hlh/hmh
5 0 1 2
t
3 hour
4
図7蓄熱量の経時変化
(粒子径の影響)
図4局所熱伝達率
Ⅷ I
b
1 1
GIQssbeods Th=12。CTt=g10。C
Prediction
一
1 −ーーー・'ロ・
5
h!f/h hlh/h
m、
mf mh
1 1
ー
1■■■
一
ィ■■■
がt=3hourでほぼ等しくなっている。これは,定常 熱伝導解より明らかにされるように,熱移動が伝導 伝熱で支配される条件では,定常状態での液相厚さ
は液相と固相の等価熱伝導率の比で決定されるため である。
対流の影響を明確にするため,平均ヌセルト数 (Nu=hmfn/スe')の時間変化の計算結果をガラス 粒子の場合に対して図9に示す。図よりDa数の増 大とともにNumhの時間に対する増加割合が増大す ることがわかる。Da=2.3×10 4 (d=16m)場合,
T=0. 1付近で準定常状態(Numf=Numh)に達し,
その後凍結量の増加,すなわち対流域の減少に伴い NumhおよびNumfは単調に減少する。d==5mの場 合,その局所熱伝達率分布は対流の特徴をよく示し ているが,その影響は小さく,粒子層内の熱移動は 伝導伝熱が支配的であることがわかる。
5. 結 言
水で飽和された多孔質層の非定常凍結熱伝達現象 におよぼす粒子径および粒子の物性の影響について 数値解析および実験により検討した結果,本研究の 範囲内で得られた主なものは次のようである。
(1) 粒子径の増大にともない自然対流の影響が表れ,
熱伝達率は増大するとともに,短時間で一定値に 到達する。
(2) 粒子層内の熱移動が伝導伝熱で支配される条件 下では,定常状態の蓄熱量は液相,固相の等価熱 伝導率比で決定されるが,粒子の熱伝導率力§大き いほど短時間で定常状態に達する。
(3) 本解析結果は,固体粒子と水の熱伝導率の比が 70以下の場合には,実験結果と良く一致する。
本研究を行うにあたり終始ご指導いただきました 北海道大学工学部福迫尚一郎教授に厚く御礼申し上 Th=12。C.1t=‑10。C,d=2mm げます。
086421QQOQ
一一Prediction
参考文献 opper
①一晩①
(1) Fukusako,S、,Seki,N.andYamada,M., JSMEht.』.,30‑264(1987),936.
(2) 稲葉,関,機論, 49‑440(昭58), 859.
(3) Weaver,J.A.andViskanta,R.,Trans.
ASME,J.HeatTransfer,108(1986),654.
(4) 岡田,冷凍, 56‑639(昭56), 3.
(5) Weaver,J.A.andViskanta,R., Int.Comm.
HeatMassTransfer,13(1986),245.
(6) Chellaian,S.andViskanta,R., Int.J.of HeatandMassTransfer,31‑2(1988),321.
(7) Brinkman,H、C、,Appl・Scient・Res,1(1947),
27.
(8) Vafai,K.andTien,C.L., Int.J・ofHeat andMassTransfer,25‑8(1980),1183.
(9) Ergun,S.,Chem.EngngProg.,48(1952),89.
mKumii,D.andSmith,J.M.,AIChEJournal, 6(1960),71.
1
UD Saitoh,T.,Trans.ASME,J.HeatTransfer, 100(1978),294.
Um Patankar,S,NumericalHeatTransferand FluidFlow,Hemisphere,NewYork.
Iron
、GI。ss
0 1 2
t
3 hour
ム
図8蓄熱量の経時変化
(粒子物性の影響)
0861 ︑0℃醐扣帥碓R
pC
o21
一一冊
DSOdqeD
SSQIG し
二Eコヱ 4︸Eコヱ
2
0 02 0A 06 08 10
て
図9平均ヌセルト数の経時変化