数理解析研究所講究録
第 812 巻 1992 年 56-66

56

体積変化と接触点の運動を考慮した
ステファン問題の数値解析
今井仁司 (筑波大学 電子・情報工学)
花田孝郎 (電気通信大学 情報工学)
河原田秀夫 (千葉大学 工学)
名取 亮 (筑波大学 電子・情報工学)

1

序

1.1
世の中の現象はすべて動的であり, その動的な現象の一つに相転移がある. 相転移は, 身近ではウィスキーのオンザロックから最先端科学においては超伝導に至るまで, 広い範囲に認められる現象である. したがって, 相転移を解析することは物理的に興味あるものであると同時に, 工学的にも重要である.
相転移現象では, 一般に, 2 つの相の間の境界が動く. 動くことのできる境界を移動境界というが, それらのなかでも境界の運動そのものが未知であるときに自由境界という. この定義に従うと, 相転移現象は自由境界を持っことになる.
自由境界を持っ現象を定式化して得られるのが自由境界問題である. 従って, 相転移現象は自由境界問題として定式される.

自由境界問題とは, 数学的には, 方程式が定義されている領域自体も方程式に依存するという問題である. たとえ方程式が線形であっても, 考えている領域も未知なので, 問題全体としては非線形問題である. そのために, 解の安定性の問題や分岐現象などが現れ, 数学的にもとても面白い問題である.
相転移にまっわる自由境界問題として, 長いあいだ数理物理学者の興味を引いてきたのが, (氷の) 融解現象をモデル化したステファン問題である.

57
1.2
ステファン問題では, 潜熱 (1 次の相転移にともなう熱) によって氷が融けるという, 熱流の平衡条件即ちステファン条件で境界の運動が規定されている[14][26]. ステファン問題の単純なものは, 空間は 1 次元で氷の領域では温度が融に設定されていて, 水の (温度が正である) 領域だけで熱伝導方程式を解けばよい場合である. これは 1 相問題と呼ばれている. っぎに, 氷解点 ($0$度とする)の部分でも熱伝れた [22]. その後, 多次元への拡張の一つに,過冷却過熱現正になることすら許容していばれ, 液相と固相の中間相と界上での表面張力) したがっ重要な因子になっている.

さらに, 平らな自由境界が不安定になる現象 [18] も見つかっている.

1.3
ステファン問題には, 以上のように様々な拡張が行われてきてはいるが, まだ満足できない点がいくっかある. 特に, 問題が非常に実用的であることを考える,と) っぎの 2 っの点がクローズアップされる(a) 相転移に伴う体積変化が考慮されていないこと,(b) 液相での流れが考慮されていないこと.

図 1: 体積変化するときの凝固現象(a) は, 液相と固相の密度差による相転移時の体積変化のことである. 冷凍庫にビール瓶を入れておくと割れてしまうといった, 水の凝固による体積膨張は日常経験することであるが, 今までのモデル化ではこのような現象は考慮されていなかった.

凝固融解の相転移における体積変化の大きな物質は少なくな58い(水で 1 割ほど) ので, そのモデル化は実用上重要であるといえる. その際に,保存しなければならない量は, 体積にかわって重さとなる. 体積変化を考慮した多次元ステファン問題のモデル化とその数値解析については, すでに行っており$[9][10][11])$ 数値計算の結果から質量保存はほぼ満足されているようである. しかしながら, このモデルにおいてもまだ満足できないところがある. それは, 容器壁の境界における流れの境界条件が都合良く設定されていることである.
緩やかな凝固的でない.

この場合には, 流れる. しかしながら, 実際には,わなくてはならない場合があ望まれる. そのためには) 流れこで問題となるのは, 固体と容点あるいは接触点と呼ばれて.

1.4
体積変化を考慮したモデルでは, 体積変化による接触点の運動の取り扱いが重要になってくる. その運動の設定の仕方によって, 領域全体における解の様子が異なってくることも考えられるからである.

図 2: 毛管現象

しかし) 接触点の問題は何もステファン問題に限ったことではない. ビーカーの中に液体を満たし, そこに細いガラス管を立てると, ガラス管の内部と外部の液体表面は一致しない (水の場合には外部よりも持ち上げられる. 図 2 参照) .
これは毛管現象と呼ばれている古典的な現象である. 現象は古典的でも, その定式化はほとんどされていない. 実際, 今までのよく知られている自由境界問題のどのモデルを適用しても, この毛管現象を動的に再現することはできない.

591.5
ただし, 静的固体の間の表面張力$\kappa$差とともにえて) 接触角$\psi$$\kappa_{W}$) を考についてのヤングの方程式$\kappa_{W}$$:=\kappa_{WG}-\kappa_{WL}=\kappa\cos\psi$が知られている [23].

静的な場合の解析はそれでよいが問題は動的な場合である. 流体は, 実際に$\{$は, 有限領域内に存在するので, 自由境界の存在や境界条件が流体の振る舞いに影響する. 従って, 実用問題として流体を解析するためには) 接触点の動的振る舞いの解析が必要不可欠になる.

2定式

なぜ今まで, 接触点の動的な解析がなされなかったのだろうと考えるに, たとえ定式化してもそれを解析する手段がなかったためだと考えられる. 実験的に,毛管内を上下する表面の様子を観察することはとても難しいし, 解析的にも困難である. 接触点の振る舞いは本質的に 2 次元以上の現象であり) しかも自由境界問題の一部であるので, 厳密に解くのはほぼ不可能であると思われる.
このように, 接触点の動的な振る舞いを定式しても, その解析は今までできなかったのであるが) 最近の数値解法の発達 [12][13][25] により数値的に可能になった.

そこで, 我々は接触点の動的な振る舞いの定式を試みようと思ったのである.

正の実数$T(\in \mathbb{R}^{+})$をとり) 領域 $\Omega=$ ]および自由境界と自由表面の規定関数ある粘性率 $1/Re$ , 潜熱形速度テンソル$\epsilon,$$\ell$$\Phi,$$0,$ $T[\cross G$

で温度$\epsilon_{ij}=\frac{1}{2}(u_{)}^{i_{j}}+u^{j_{i}})$, 流速$u$, 圧力 p)を未知として考察する. 物理定数で$\Phi^{+}$, 比熱 , 熱伝導率 , 表面張力$c$$\theta$$k$$\kappa$

および平均曲率 $H$ , 変とともに方程式系は無次元化して表現する.

60

2.1る固体領域

温度は氷の$\Omega^{s}$

では熱拡散方程式

in$c\dot{\theta}+k\nabla^{*}\nabla\theta=0$を) 温度が正である液体領域$\Omega^{L}$(1)$\Omega^{s}$

では流速を含んだ拡散方程式$c \dot{\theta}+k\nabla^{*}\nabla\theta-c\nabla^{*}(u\theta)-\frac{1}{Re}(\nabla u)\cdot\epsilon=0$

in$\Omega^{L}$(2)を満足している. ここで記号は時間微分,間微分
gradient,

い篭$-\nabla^{*}$

は発散 divergenceである. 液相は非圧縮性の粘性流体と考えるので連続の式

in$\nabla^{*}u=0$(3)$\Omega^{L}$と Navier-Stokes 方程式$\dot{u}+(u\cdot\nabla)u=-\nabla p-\frac{1}{Re}\nabla^{*}\nabla u$in$\Omega^{L}$(4)を満足する.

2.2の境界上で満たすべき条件を列挙する. 容器は底からっぎに, 領域だけ冷却されていて, 壁では断熱的と考えれば$\Omega,$$\Omega^{S},$$\Omega^{L}$$\theta$$=$$\theta_{\Gamma}(<0)$$\theta_{\nu}$

$=$

$0$

on

$\Gamma^{W}$

on ,$\Gamma$

(5)

(6)

61
を, 液相での流れは壁に沿ってすべる条件$\nu\cdot u$

$=$

$0$

$(u_{\mathcal{T}})_{\nu}$

$=$

$0$

,(7)

on

(8)

$\Gamma^{L}$を仮定する.
液相と固相の内部界面\Sigma$=\partial\Omega^{S}\cap\partial\Omega^{L}$

上では温度の連続性(9)および Stefan 条

(10)
および

(11)

$(1-\alpha)\dot{\Phi}=-u\cdot\nabla\Phi$を仮定する. ここで定数$\theta^{L}=\theta|_{\Omega^{L}},$$\alpha=\rho^{S}/\rho^{L}$$\theta^{s}=\theta|_{\Omega^{S}}$,は膨張係数 (密度比) であり,$k^{L}=k|_{\Omega^{L}},$$k^{s}=k|_{\Omega^{S}}$

,

である.

自由表面\Sigma +$=\partial\Omega^{L}\backslash (\Gamma^{L}\cup\Sigma)$

上では

(12)

$\theta_{\nu=0}$を仮定する. また, 自由表面の保存条件および応力の連続性から$\dot{\Phi}^{+}$$\nu\cdot\epsilon\tau$

$p$

$=$

$-u\cdot\nabla\Phi^{+}$

,

(13)

,(14)

$=$

$0$

$=$

$2 \kappa H+\frac{1}{Re}\nu\cdot\epsilon\nu$(15)を満足する.

最後に, 接触点付近での自由境界における力のつりあいを考えて, 接触角に関連する条件を導入する [1]. そのときには, 境界層のようなものの存在を仮定しておくと$u^{z}=Re(\kappa_{W}-\kappa\cos\psi)\tan\psi$on$\partial\Sigma^{+}$(16)62が得られる.

速度の境界条件として粘着条件を用いるには, 整合条件 (compatibility condition) を仮定するかどうかで状況を分けて考える必要がある. 整合条件を仮定する場合は, 粘着条件を満たしている境界部分と接触点の間に) 速度を滑らかにつなぐ部分境界を考える必要がある. これは遷移境界 (領域) と呼ぶが, そこでは例えばスリップ条件を仮定することになる. 整合条件を仮定しない場合には, 特に付加条件を考える必要はない

3数値計
数値的には)]. 空間の離散化は差分法 (移(速度場温度場と写像関数を交互に求める方法) で行う.

物理定数としては $k^{s}=24,$ $c^{s}=1.2,$ $k^{L}=5.6,$ $c^{L}=5.6,$ $Re=20,$ $\alpha=$

0.917, $\ell=100,$ $\kappa=100$ として, さらに物理的な考察 ( $\cos[\psi_{0}]=0.2$ なる平衡接触
角 $\psi_{0}<\pi/2$ として) から式 (16) は$u^{z}$$=$$Re(\kappa_{W}-\kappa\cos\psi)\Psi[\psi]$$\Psi[\psi]=\{\begin{array}{l}tan[\psi])tan[\lrcorner\end{array}$on$\partial\Sigma^{+}$,$if\psi\leq(.<if\psi>^{\frac{\psi}{\underline\psi}}\frac{\pi}{2})$初期状態で水と氷の深さと半径の比は 0.2:1:1 で, 底と上表面での温度の比を $-2:0.01$ として (図 5) 計算を始めると図 6 の状態まで結果が得られてとしている.

いる.

参考文献
[1] $G.I\langle.Batchelor$ : An introduction Fluid Dynamics, Cambridge Univ. Press,
1967[2] M. Bertsch, P. de Mottoni, L.A. Peletier: Degenerate diffusion and theStefan problem, Nonlinear Anal., 8(1984), 1311-1336.
[3] M. Bertsch, P. de Mottoni, L.A. Peletier: The Stefan problem with heating;Appearance and disappearance of a mushy region, Trans. Amer. Math. Soc.,
239(1986), 677-691.

63
Tlme
otal$T$

$=$

mass
$=$

Ice

1.102294$regi=^{on}0.9$

Water

02041
0200253

region$=$

図 5: 初期状態

0.000000

64TimeTotal

$=$

mass
$=$

Ice region
$=$

Water

図 6: 中間状態 $(t=0.04)$

0.040000
1. 102363
O. 919053

region$=$$01$ 83309

65[4] G.Caginalp: An analysis of a phase field model of a free boundary, Arch.
Rational Mech. Anal., 92(1986), 205-245.

[5] J.R.Cannon, E.DiBenedetto: On the existence of weak-solutions to an ndimensional Stefan problem with nonlinear boundary conditions, SIAM J.
Math. Anal., 11(1980), 632-645.
[6] J.Crank: FrdibdblCldPess, Oxford,1984[7] A.Fasano,ions of onedimensional), 115(1977),

341-348[8] A.Friedman: Te Stefan proem in severa variables, Tlrans. Amer. Math.
Soc., 133(1968), 51-87[9] T.Hanada, H.Imai, H.Kawarada, M.Natori: Numerical Computations forSolidification problems with Change of Volume, Bull.of the Greek Mathe.
Soc.,31 (1990), 29-49[10] 花田孝郎, 今井仁司, 河原田秀夫, 名取亮: 凝固現象における自由境界問題の数値解析, 数理解析研究所講究録, 744(1991), 100-111[11] 花田孝郎, 今井仁司, 河原田秀夫, 名取亮: 体積変化を伴う凝固現象の数値解
析, 数理解析研究所講究録, 746(1991), 116-129[12] F.H.Harlow, J.E.Welch: numerical calculation of time dependent viscousincompressible flow of fluid with free surface, Phys. Fluids, 8(1965), 21822189[13] Y.Katano, T.Kawamura, H.Takami: Numerical study of drop formationfrom a capillary jet using a general coordinate system, Theor. Appl. Mech.
34(1986), 3-14

[14] H.Kawarada: 自由境界問題, 東京大学出版会, 1989
[15] 村田健郎, 名取亮, 唐木幸比古: 大塑数値シミュレーショ ン, 岩波書店, 1990
[16] M.Natori, H.Kawarada: 自由境界問題の数値解析, Japan. Phy. 31(1976),
547-551

66[17] R.Kobayashi: A mathematical one-dimensional one-phase model of supercooling solidification, Publ. Res. Inst. Math. Sci., 19(1983), 327-344[18] C.W.Lan, S.Kou: Thermocapillary flow and natural convection in a melt
column with an unknown $melt/solid$ interface, Int. J. for Numer. Met. inFluids, 12(1991), 59-80
problems us[19] R.H.Nochet
ing numeri84-814[20] R.H.No chetate parabolic
Math.Comp.
free bounda
51(1988), 2[21] T.Nogi: A methematical one-dimensional model of supercooling solidification, Publ. Res. Inst. Math. Sci., 21(1985), 1121-1203[22] O.A.Oleinik: A method of solution of the general Stefan problem, Sov. Math.
$Dold.1)$ (1960) 1350-1354[23] 小野周: 表面張力, 共立出版, 1980[24] D.Takahashi, Y.Takeda, H.Takami: Numerical simulation of collision ofliquid droplets, Theor. Appl. Mech. 36(1988), 3-15
[25] Joe F.Thompson, Z.U.A.Warsi, C.Wayne Mastin: Numerical grid generation, North-Holland, 1987

[26] 山口昌哉, 野木達夫: ステファン問題, 産業図書, 1977