NAOSITE: Nagasaki University's Academic Output SITE Title デルタ列関数によるステファン問題の解析 Author(s) 福地, 信義 Citation 長崎大学工学部研究報告, 15(24), pp.59-64; 1985 Issue Date 1985-01 URL http://hdl.handle.net/10069/24177 Right This document is downloaded at: 2014-11-01T13:20:11Z http://naosite.lb.nagasaki-u.ac.jp 長崎大学工学部研究報告第15巻第24号．昭和60年1月 59 デルタ列関数によるステファン問題の解析福地信義＊ Analysis of Stefan Problem by Using the Function of Delta−Sequence by Nobuyoshi FUKUCHI＊ The oil tankers￠arrying high pour point crude oil almost are equiped with the tank heating system consisted of meandering steam pipe to reduce the．viscosity of crude oil． However， the accident of oil solidification occasionally occurs during pumping up cargo o量l from an oil tank of tanker to land for the sake of．unsuitable oil heating system． For a settlement of this case， the heat energy of steam is transfered to crude oil in solid state to melt it． In．ord6r to量nvestigate the heat transfer with moving．solid phase boundary， so・called Stefan problem，・ two dimensional numerical analysis is carried out in which the enthalpy is expressed continuously by us童ng the function of delta＿sequence at interface of both phases． Furthermore， the spreaded melting zone and the oil temperature by using paraffin are measured in quasi−two dimensional tank with two heat3rs． The results of theory and experiment are compared for judging the lustice of analytical method． 1．はじ．めに解析し，固液両相内の温度分布と液相内の対流速を計高流動点原油を輸送するタンカー．において，その装算する。また凝似2次元タンクによるパラフィンの溶備した加熱システムが適当でない場合には船底部に原融実験を行い，溶融域の進展状態，温度分布などにつ油を凝固させ，揚油不能となる事故が起ることがある．いて数値計算の結果との比較を行い，解析法の妥当性この場合には，船底付近に設置された蒸気管に通気を調べる．することにより凝固状態の高流動点油を溶融し，さらに低粘度化を行なう． 2．ステファン問題このような，固液相境界が移動する場合（わ熱拡散， 2．．1 基礎式いわゆるステファン問題1）を解析するためには相境界高流動点油を加熱により溶融する場合のように，相における融解潜熱の扱いが問題となる．．本論文では相変化があり，固液相境界が移動する場合の熱拡散の状変化時におけるエンタルピーをデルタ列関数を用いて態は以下の諸方程式により表わされる．連続的に扱い，固定メッシュの有限要素法により数値（1）液相内昭和59年10月1日受理＊構造工学科（Department of Structural Eng量neering） 60 デルタ列関数によるステファン問題の解析液相内では自然対流による熱伝達が起こる．油を近 2．2 デルタ列関数を用いた熱拡散方程式似的にニュートン流体とし，流れを層流とすると，速（1）エンタルピー関数度場と温度場の状態を表わす式は次の通りである．（5）式で表わされるステファン条件を相境界で厳密なお式はtensor表示とし，添字は総和規約に従うものに満足させ，熱拡散方程式を解くことは困難なので，とする．次のエンタルピー関数H（θ）を用いて近似的に解析を 1）流れの連続式行う2 H（θ）イ・（θ）μθ）dθ＋μθ・）仇Y（θ一θ・）（6）器一・（δ＝1，2，3）（1）ここに，π‘は流速である．oσ妥は座標を表わし，」σ、，コσ2 ここに，Y（コσ）はHeavisideのステップ関数であり，を水平方向，」C、を鉛直方向（上向きを正）とする． θmは油の融点である。 iD運動方程式高流動点原油のような混合体では融点にある程度の重力場にお・けるNavier−Stokesの運動方程式は次式幅があり，また比熱，比重量の温度依存性が小さいことで表わされる．を考慮するとエンタルピー関数は次のように表わされ讐＋毒（・幽）一一去藷る。＋素｛・（∂％‘十∂賜∂：r」 ∂：r‘）｝一β9・θ（2） H（θ）一・ρθ＋呵コω・（θ一θ・）dθ（7）ここに，Pは圧力，孟は時間，θは油温である．ここに，ωk＠）はデルタ列関数であり，Diracのデまたg、＝g2＝0， g3＝一g（重力の加速度）であり，レルタ関数δ（勾との間には次の関係がある．は動粘性係数，βは体積膨脹率，ρは比重量である． δ（oσ）＝limωk（oじ）（8）（2）式の境界条件としでは固液相境界およびタンク周 k→oσ 壁においてno−slipの条件をとり，自由表面では鉛直デルタ列関数の例を3種類Fig．1に示す。ここでは方向の流速を0とする・最も扱い易い， liD熱拡散方程式 ωk（oo）＝κ／π（1十κ『コσ2）（9）器＋素（ω，θ）「妾（畷）（3）を用い，κとωk（00）との関係およびこれに対応するエここに，κ、（＝λ、／cρ）は液相内の熱拡散係数，λ・ンタルピー関数正1（のをFig．2に示す。デルタ列関数をは熱伝導率，cは比熱である．用いることにより，融点において潜熱のため不連続と（2）固相内なるエンタルピーを連続的に扱うことができる．凝固相の熱伝導を表わす式は（3）式において移流項がない次の方程式である．誓一素（ ∂θκ8 ∂：r」）（4） Wk 亜 30 L一冷㎏蝉 2．5 ここに，κ・（＝λ、／cρ）は凝固相内の熱拡散係数である． 20 （3）固液相境界相境界では，流出入する一流束の差が融解潜熱に等 t5 しくなる，いわゆるステファン条件と呼ばれる次の熱 1．0 ’ 収支方程式を満足しなければならない・擁論、レ響 ∼ ’ 為（∂θ∂η）。一為（誓）、一ρσ・誓（5） 1 Q5 、1 1’ ここに，ηは固相境界に対し外向き法線方向を表わ 0 積当りの融解量を意味する．またσ、は融解潜熱である．一3 、・乏 ∀1 ）し，添字8，Zをそれぞれ固相側，液相側の値を示す．呪は・η方向の境界移動量，すなわち境界上の単位面へ／へ一〇，5 Fig．1 ‘’ 、 ’・〆・、 1 》冨・3 0 Three examples of function of delta− sequence。福地信義 35 60 Wk 誰Wk＝禰 30 （x104） r票∫／詮50 ll／／ H（θ） ii ii 1i HeQV 54ρ ／ oil @ イ口 O0 （200 む ii 20 RユrQ箔 O0 妻45 lレk＝10 2．5 00 ．鋼艦：翌一 55 61 緬ノ 1．う ii 35 i：k．5 ao ／H・c・θ＋q・践 t5 25 20 1．0 k＝2 Q5 ゆツ 15 k＝1 ／1 ．／1 ；100 O0 3 ヤムへ暮／亀．逃〉ぐ T0 ’ぬ。 lo目e oil 勤亀＼ Mi @ ノ＼ pS ／／ソ1 11 ＿7豪ニーメ’ @ ，ノノ Q0 モ＼ р・@oi｛、℃ @ で）［ず一＿ tO 一一一 0 受心、，’窪，ノ melting point O5 一4 −2 0 2 4 x 3D 10 ﾍt轟，．，罰θおec） Fig．2 A fUnction of delta・sequence and enthalpy 噂 ShQi ・くy ，Rユ十N十†十手「＼辱 ffin B貢モ窒浮 @5 ゆ十 ● 十十幸「一、 fUnCt孟On USing it 210 2Q 30 40 50 60 70 80 temperQture of oil（。C）〔fil㎏d up by pαrqffin〕 ceUing：（1cryb乏e …m／ †，・。、d認聰） item uhit symbo1 丁ach拍g モ窒浮р・@oi1 Minas モ窒浮р・@oi1 1 （8柵kne＄）暢〕層ゴ漏 rqd iQtive 卜edねr セGnslbrmer thermoイbw meter チ bottom：「x）1ystyro【 fOQm m・Qsu…g・㎝・・・…．門ヒ…一r T initial using sur口 heGter 秩Bunding 盾奄s ・・≠狽・ code （34．2W）狽・撃o．狽・香D 陛議 P。ur p。int o C Melting point 3Kg／m Paraffin 浮唐・п@for ・・@erlment o C θP 30 35 鼻一一臼5 30．86×10 30．85×10 一 42’》44 30。81・X10 Speciffc weight 3Kg／m ρ ρ・9，（1一・．7・1・一3日9 ＝ 15 （1−0．9Xτ09 Cubical ・・垂≠獅唐奄盾氏@rate Specific heat 1／K 」／Kg・K Latent heat 盾・@fusion Themal モ盾獅р浮モ狽lty ・」／K ．W／m・K β C qL λ 0．57／ρ C司ア40／／繭δ＋3．8（θ一15） 2．。9∼2．26×105 λ・117（1一．。．54・10−3θ）／ρ ・・t・・θ・．・i1・・mper・t・re（．C）・θ15・θ一15 （right）35．9W M2−1 28。7。C 21．3。C i1eft）34。2W M2−2 36．1 26．0 M2噂3 42．5 19．6 Fig．4 Physical properties of high・pour piont o三1 （unit：mm） Fig．3 Quasi−two dimensional tank model and ここに，D／D孟は流体力学的微分であり，1固相では measuring condition． D／1冗＝∂／∂ち液相ではD／D置＝∂／餌十篤∂／∂筋となる．さらに（11）式は次のように書き改められる．なお，（9）式におけるκの値｝よH（θm）＝ρθm（C＋αL／［1＋σLω。（θ一θ。 C）P器一蹴（・署）（12） θm）とωhの最大値ωk（o）ニκ／πとの比較からκ＝πσL／ θmとした．（2）熱拡散方程式 3．解析結果（7）式のエンタルピー関数を顕熱容量で割ると次の 3．1 パラフィンによる溶融実験2｝ようになる．（1）実験方法．一ん（θ）一 ∴齠ｻπω・（θ一㊧dθ（1① 高流動点油の加熱による溶融範囲の拡がり方と温度分布を調べ，．さらに2．．2で述べたデルタ列関数を用このん（θ）を用いることにより，、（3）式および（4）式で表わされる熱拡散方程式は，相境界における潜熱吸収の取いた解析法による解とを比較するために，ガラス製のし凝似2次元タ．ンクに融点42∼44℃の固形パラフィンを扱いを考慮した，すなわち弱い意味においてステファ入れて，2本の円柱状電気ヒーター（加熱管と称する）ン条件を満足する次の方程式に変換される．により加熱を行い計測を行った．この実験装置と実験号｝一素（ ∂θκ ∂lr」）r．（11）要目をFig．3に示す．溶融域の拡がり方はガラス面に貼ったフィルム上にデルタ列関数によるステファン問題の解析 62 相境界を直接記録した。温度分布は非接触型熱放射計考えられる．を用いて，縦横11×13の格子点上の熱放射量を計測し，速度と温度に関する方程式（1），（2），（12に対しGaler・経験的に求めた熱放射量一温度の関係式により温度へ kin法を適用して有限要素法のための定式化を行い，換算した．なお熱放射量は周囲の熱的環境に影響され 2次元問題として扱った．有限要素としては8節点のるので，これを考慮して温度への換算式を作った． 2次要素を用い，定式化後の矯，P；θに関する非線形実験に用いたパラフィンの粘度一温度特性およびそ連立方程式はNewton−Raphson法により数値的に解き，の他の物性値を高流動点原油の代表例であるミナス原時間微分項は差分で近似した．油（インドネシア産）および大慶原油（中国産）と比計算対象はミナス原油とし，Fig．7に示す計算条件較して3）・4）Fig．4に示す。により解析を行った．また加熱開始時には加熱管まわ（2）溶融状態と温度分布りの3×3の要素は液相要素とし，他を固相要素とし加熱時の溶融域の進展の様子をFig．5（a）に示す．た．固相要素はその平均温度が流動点以下であるか否なお，Fig．5（b）は，3．2で述べる数値解析法によりかを判定し，超えている要素は液相要素に転換した．求めた溶融域の時間的変化である．溶融は最初加熱管（2）溶融状態・上方に向って起こり，上面に到達すると液相が固相よ液相内対流について計算した流速分布の時間的変化り軽いことと膨脹により上面横方向に液相が拡がり，の様子をFig．7に示す．対流は加熱管上方への上昇流次に上面から下方に向って溶ける．加熱管より下部はと側方の下降流から成り，溶融域の広さに応じて循環対流が起らず移流の効果がないため，伝熱量が少ない流の数が変化するのが見られる．ので，加熱管より上部が完全に溶解するまで融解しな Fig．8はF量g．7に対する温度等高線の時間的変化である．パラフィンによる計測でも見られるように，固い． Fig，6は弾弓分布の経時変化を示す．いずれも固液相境界では融解潜熱が大きいために，両相とも境界に対し法線方向の温度勾配が大きくなるのが見られる， 3．2 数値解析（1）計算法相境界が移動する場合の熱拡散状態を解析する方法として，縦横8×9の固定メッシュによる有限要素法を用いる．ただし，固定メッシュのため移動する境界 ②竜一120mln を完全には追従できないが，相境界付近では非常な高 ⑤t・330min， 40 粘度のため流れはほとんど起らず，流速計算において弓。境界位置のわずかな食い違いあまり問題にならないと 60 tOP 60 50 6 4 む 15018・ 45 1：。 30 4 ；：1 ③t−180mln 90、 15 15 1i。 ⑥t・420min．：1： 8 270 1：よ 27 @8 0 30010 ll。 330 00 ii§ P2 iii 330 il崔 420 heαter meGsurement M2＿t ㎞甘 12 bottom （unit of time ：minute）（G） π灘3urernent h・・tq・antity・Q1（θH一θ1階12x2 （・nit・f ti…λ・／・戸2・。．、・・1・一3）（b） CGIcutGtion （unit：。C ） Fig．5 Growth of melting zone in tank model Fig．6 Distribut量ons of oil temperature in with two heaters． measurement M2−1 福地信義 63 ρ。0．ヨx1σ3 ④入tlc ①撞ノ。 28 9 【2。2鴉x1σ3 0 3 3 3 8 4 37 0 38 4 一 36 入t1・ρ匹147・1σ3 ⑤入ゼcρ匹331x1σ3 7 2 2 7 8 ② 入t1、12．11。．1σ3 ⑤入t1、12。3．31．1♂ 4 52 2 ロロロロ 54 44 4 34 ・L… 亀■←隔。79 噛し■ … ．．・艫潤@コべいやぐなゆかルジロロコヒ：；1獄でな駕ニモ㌧二＝；：： i癩鱗髪磁；ii：・・ 46 56 6 ロののうちもロキロヲロもメのロほコ 40 3 3 58 0 ’曹召’ @ゆゼ民」rガ． ○、＼庵’メ」・や・・ ’・ @・腎¶ ．，o・9摂 r r も煽○貝ρ，コやロ 3 5 ゆサロやひコそセゆロサリヘヨぱサロロげヒ ’：；＝や、こ＝：：：：：1二二三ナこ：：：ココロものぼロコロヒコロリロロロサヒヒコ：：：：∵：：：：：：：：：：∵｝：：，： 8 ③入tlcρ【≧220x1σ Calculating condition of ta6k mode1 2 3獄，c l2．1．脳．1（戸（unit・・C） 1） Assumed length ：1旨 2 閉 heGt quqntity 2） Heating oi1 ＝Minas crude oi1 （pour point ： 340C ） @ Q／（θH一θ1）KA＝6x2 3）Oil temperature near heati・g・・i1・θH・6・．C 8 4）011t・mper珈re at b・gimi・g・θ1・3・．C 40 日ote…． j：Heat・transm1ゴsi。n c。efficient 5）S・一di・g t・mper・ture・醜・25．C 36 ‘）日eat B辮翻llgC3：｛1’鷲．8、，．2K @A ＝ Heat transmission ared i…τ・mper・t・re・e・r heati・g・・i1 ﾆ1・Oil t・mper・加爬・t b・gimi・g 下 6 0．6875 275 Q275 掴45奇 Q3125 ． heαter ⊥ Fig．7 Calculated flow pattgrn in model C2−2 Fig．8 Distributions of calculated oil temper・ and calculating condition． ature in model C2＿2 1．0 qL冨0 qL＝218（J／㎏）弩お04 一L． latent heat M2〆 total 窪1・0 9・6 M2−2（36．10C） 9Q8 M2−3（42．50C） 9 M2−1（28．70C）歪oβ te「nperoture 、 surrounding ち04 220．2 ε he（1t qLkユntity 2Q2 Q／（θH一θ【）KA＝6x2 0 （G）meQSUrement qL＝1・09 U！㎏）轟・・ 0 ．Q5 10 C、㎜t ｯ。轟、225 01234567 0 3＆σ・）蕊 8 9 10 セime（hour） Fig．9 （b） CGヒu［otion Variations of calculated melting zone 1．0 岩、， corresponding to three kinds of latent heat． 9 婁Q6 ＼も04 液相境界において両相の温度勾配は融解潜熱吸収のため大きい．これに反し，境界から離れた固相部では伝熱量が少ないため，また液相部では対流による掩乱のため温度が均一化し，温度勾配が小さくなっている．＼1；二1羅｝iW ￡ 2Q2 0 C2−1（q＝3x2）ロ note二q＝Q！織一θp KA 0 α51 唐謔T．．遣2，、225 30 35 （x1σ3） 2．2で述べた，デルタ列関数を用いた熱拡散方程 Fig．10 Variations of melting zone：correspond・式により相境界の熱収支において潜熱の扱いが妥当で ing to supplied heat quantity．デルタ列関数によるステファン問題の解析 64 あるかを調べるために，高流動点原油の融解潜熱σL＝ 4．結言 2．18×105（J／kg）に対しその％および0の場合について解高流動点油を対象に；凝固状態から加熱により溶融析を行い，この3種類の場合についてタンク全体を1．0 する際の相変化を含む熱拡散の状態を解析するために，とした溶融量比の経時変化をFig．9に示す．溶融量比熱拡散方程式においてエンタルピーをデルタ列関数をの曲線は固定メッシュのために波うっているが，潜熱用いて温度に関する連続関数として表わし，固定メッの大きさの差による溶融速度の違いが現われている．シュの有限要素法による数値解析を行った．またζれ Fig．10（a）はパラフィン溶融実験にお・ける外周温度に対応するパラフィンの溶融実験を行って解析解とのと溶融量比との関係であり，加熱開始2∼3時間は初比較を行った．これにより完全ではないが本解析法に期油温に影響されるが，その後は周囲に逃げる熱量のよるステファン問題の解析が可能であることが分った．多少により溶融速度が決っている．これに対して，外しかし，溶融域の進展が階段状になるなど，移動境部への貫流でなく供給熱量を変化させた場合について界の取扱いやデルタ列関数におけるκ値の選択法など計算した溶融量比の時間的変化をFigユ0（b）に示す．に問題が残されている．固相と液相およびタンク全体の各平均温度の経時変本研究は文部省科学研究補助金を得て行ったことを化を，計測値と計算値を比較して，Fig．11に示す．タ付記する．ンク全体の平均温度は溶融域が小さい初期には固相の温度に近く，溶融が進むに従って液相の温度に近づくが，加熱管より上部が全て溶け溶融速度が鈍ると潜熱参考文献に奪われていた熱が液相温度の上昇に使われ，油温が 1）山口昌哉，野木達夫：ステファン問題（1977），急激に上昇する．固相温度の時間的変化は小さく，こ産業図書の図からも伝導のみによる伝熱は少ないことが分かる． 2）福地信義：高流動点油の加熱による溶融と熱拡散について，日本造船学会論文集，第156号（1984）（1） meqSUrement 3）渡辺，他：貨油タンクの高粘度原油洗浄，昭58春 60 70 季船舶技術研究所発表会講演集（1983）鯉P蟹ノ 55 65 ノト虚』’ルぜ 0 劉量45 4）日本造船研究協会：タンカーのタンクヒーティン月目。／＿ノ説・鏡謹議／ヨ量55 塁、。警〆崔go Aノ簿蓉35 45 一一鰻鰹∼’ 30 グに関する研究，協会資料，No．79（1968），No．99（1969）ゲ／び 40 25 35 20 30 15 25 99［id phqse ／阜一｝響，・・ひ●一一 01234567 0 1 2 3 4 5 time（採）url （b）meosuremenい42−2 tr爬（㎞r）（Q）meqsurement M 2−1 （2）co【culotion 48 50 48 Heat quantity Gy（像一θ，）漁＝6x2 ε46 46 44 睾44 42 §42 540 ε 1iquid phase 豊暴2 譜ｪび 36 ’豊 Heat quantity Q／（一一a）KA＝12x2 評5e 、、Φ、も 40 二 538 ，〆 38 説漣 36 巨 9 34 34 ＝ 32 SOlid phase o 32 28 solid phase 30 30 0 0．5 1。0 1．5 @ time ．二又，，、誕5・・1δ9・・ 28 むらしむの段1げ・｝2’‘ 、b、、。。，u、謝謝（o） ooヒu1αtion l C2−1 Fig．11 Variations of averaged oil temperatures in solid phase and liquid phase．