自由表面流れにおける代数的界面再構成法であるTHINC法を紹介し，その計算過程で現れる不適切な結果を，物理的観察で回避する方法を述べます．
自由表面流れの計算法
流体とは，気体と液体の総称です．液体がある種の容器に貯められているときに，液体とその上方に存在する気体の界面を，自由表面とよびます1．海面は我々がよく知る自由表面の代表例です．
この自由表面を伴う流れ（自由表面流れ）の数値計算は，船体抵抗の予測，津波の遡上などの観点から多く行われています．自由表面流れの数値計算法はいくつかに分類することができ，その中でもよく使われているのが，1流体モデルと体積率に基づく計算手法です．
1流体モデルでは，気体と液体を1本の式で表し，液体の場所・気体の場所で物性値（密度および粘度）を切り替えます．液体・気体の識別や，流れの時間発展に伴う液体・気体の場所の更新には，Volume-of-Fluid法（略してVoF法）を用います．VoF法では，VoF値とよばれる液体の体積率（液体の場所で1，気体の場所で0になる関数の局所的な空間平均）の空間分布を求め，その空間分布を次の移流方程式に基づいて移動させます．
\frac{\partial \phi}{\partial t} + \nabla\cdot(u\phi)-\phi\nabla\cdot u=0
$\phi$は液体の体積率，$u$は速度です．
VoF法の考え方はとても簡単でわかりやすいのですが，難点は，移動後のVoF値からの自由表面の復元です．ある微小な体積（セル）の中でVoF値が0.7と分かったとして，自由表面がどのように分布しているかを推定するのは簡単ではありません．ある体積内のVoF値だけでは推定できないので，例えば，周囲のVoF値も用いて自由表面の法線ベクトルを計算し，表面の形状を推定します．このように，VoF値から自由表面の形状を推定する過程を，界面再構成とよびます．以降，自由表面を界面と称します．
VoF法は，界面再構成の種類によって，幾何的VoFと代数的VoFに分類できます．代数的VoFは，界面の形状を関数（界面再構成関数）で表現する方法です．代表的な方法に，THINC(Tangent of Hyperbola for INterface Capturing)2があります．
THINC法では，その名前の通り，VoF値が0より大きく1より小さいセルの中において，界面の形状を$\tanh$で表します．
1次元で説明すると，流体体積率の移流方程式
\frac{\partial \phi}{\partial t} + \frac{\partial u\phi}{\partial x}-\phi\frac{\partial u}{\partial x}=0
において，セル$i$におけるVoF値を，次のように$\phi$の空間積分として表します．
\Phi^n_i=\frac{1}{\Delta x_i}\int_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}\phi(x,t)dx
$\Phi= 1$のとき，そのセルは液体であり，$\Phi= 0$のときは気体を意味します．$0<\Phi<1$のとき，そのセルに気液の界面が存在する事になります．
THINC法では，界面再構成関数を次のように定義します．
f(x) = \frac{1}{2}\left\{ 1+ \alpha \tanh \left(\beta\left[X-d_i\right]\right)\right\}
ここで，$X$はセル（$x\in[x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}]$）内で規格化された距離（$=\frac{x-x_{i-\frac{1}{2}}}{\Delta x_i}$），$d_i$は界面位置（実際には$f=0.5$となる位置）です．$\alpha$は界面の向きによって変化するパラメータであり，$\phi_{i-1}<\phi_{i+1}$のときは$1$，$\phi_{i-1}>\phi_{i+1}$のときは$-1$をとります．1次元軸に対して垂直な界面の（液体向きの）法線ベクトルのようなものと考えられます．$\beta$は，界面の厚みを調整するパラメータです．
つまり，セル内において，体積率$\phi$の分布を$f(x)$で近似するので，VoF値は
\Phi^n_i=\frac{1}{\Delta x_i}\int_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}f(x)dx
として計算されます．
THINC法では，移流項の計算をセルの両界面を通過するFlux $F$の差として計算します．
\Phi^{n+1}_i = \Phi^n_i-\frac{F_{i+\frac{1}{2}}-F_{i-\frac{1}{2}}}{\Delta x_i}+\Phi^n_i\frac{u_{i-\frac{1}{2}}-u_{i+\frac{1}{2}}}{\Delta x_i}\Delta t
Flux $F$は，例えば$x_{i+\frac{1}{2}}$における速度$u_{i+\frac{1}{2}}$が正のとき，次式で計算します．
F_{i+\frac{1}{2}} = \int_{x_{i+\frac{1}{2}}-u_{i+\frac{1}{2}}\Delta t}^{x_{i+\frac{1}{2}}}f(x)dx
時刻$t$におけるVoF値が分かっているとき，THINC法は次の手順で時刻$t+\Delta t$のVoF値を得ます．
- VoF値$\Phi^n_i$と近似関数$f$の関係から$d_i$を算出する．
- Flux $F_{i+\frac{1}{2}}$を計算する．
- VoF値の移流方程式を計算し，$\Phi^{n+1}_i$を得る．
不適切な結果
THINC法では，手順1の界面位置を計算する際に，不適切な結果が得られる場合があります．THINC法の計算過程を見たとき，「気体あるいは液体のセルにおいて，界面位置$d_i$はどのように表されるのだろう？」と思ったことでしょう．不適切さは，この疑問と関係します．
まず，界面位置$d_i$を求めていきましょう．簡単のために，$x$を規格化された$X\in[0,1]$で表現し，$\alpha=1$とします．このとき，VoF値$\Phi$は次式で表されます．
\Phi_i=\int_{0}^{1}\frac{1}{2}\left\{ 1+ \tanh \left(\beta\left[X-d_i\right]\right)\right\}dX
これを積分すると
\begin{eqnarray}
\Phi_i&=&\left[\frac{1}{2}\left\{ X+ \frac{1}{\beta}\log\left( \cosh \left(\beta\left[X-d_i\right]\right)\right)\right\}\right]_0^1\\ &=& \frac{1}{2}\left\{ 1+\frac{1}{\beta}\log\left( \frac{\cosh(\beta-\beta d_i)}{\cosh(\beta d_i)} \right) \right\}
\end{eqnarray}
なので，
\beta[2\Phi_i-1]= \log\left( \frac{\cosh(\beta-\beta d_i)}{\cosh(\beta d_i)} \right)
から$d_i$を求めることができます．
両辺の指数をとり，$\cosh$と$\exp$の関係式に従って変形します．
\begin{eqnarray}
\exp\left(\beta[2\Phi_i-1]\right)&=& \left( \frac{\cosh(\beta-\beta d_i)}{\cosh(\beta d_i)} \right)\\&=&\frac{\exp(-\beta+\beta d_i)+\exp(\beta-\beta d_i)}{\exp(-\beta d_i)+\exp(\beta d_i)}
\end{eqnarray}
ここで，記述を簡略化するために，次の変数を導入します．
\begin{eqnarray}
e_1 &:=& \frac{\exp(-\beta+\beta d_i)+\exp(\beta-\beta d_i)}{\exp(-\beta d_i)+\exp(\beta d_i)}\left(=\exp\left(\beta[2\Phi_i-1]\right)\right)\\
e_2 &:=& \exp(\beta d_i)\\
e_3 &:=& \exp(\beta)
\end{eqnarray}
$e_1$を整理し，$e_2=$の形に変形します．
e_1 = \frac{e_2/e_3+e_3/e_2}{1/e_2+e_2} = \frac{e_2^2+e_3^2}{e_3(e_2^2+1)}
e_2^2 = \frac{e_3(e_3-e_1)}{e_3e_1-1} = \exp(2\beta d_i)
両辺の対数をとり，整理することで，$d_i$が得られます．
d_i = \frac{1}{2\beta}\log\left(\frac{e_3(e_3-e_1)}{e_3e_1-1}\right)
Φ=1のとき
$e_1 =\exp(\beta)$となるので，$e_3-e_1=0$から，$d_i=\frac{1}{2\beta}\log(0)\rightarrow-\infty$となります．
Φ=0のとき
$e_1 =\exp(-\beta)$となるので，$e_3e_1-1=0$から，0除算が生じます．
よくある対策
疑問に思った通り，気体あるいは液体のセルにおいて，界面位置$d_i$が$\infty$あるいは$-\infty$となり，結果としては不適切です．この対策としてよく行われるのは，微小な値$\varepsilon$を加え，$0$除算や$\log(0)$を回避することです．
d_i = \frac{1}{2\beta}\log\left(\frac{e_3(e_3-e_1)}{e_3e_1-1+\varepsilon}+\varepsilon\right)
しかし，$\varepsilon$にどのような値を用いればよいのかは不明であり，計算が発散しない値を試行錯誤的に探っているのが実情です．また，$\Phi_i=0$のとき，$d_i=\frac{1}{2\beta}\log\left(\frac{e_3(e_3-e_1)}{\varepsilon}\right)\rightarrow\infty$となり，界面位置として$\infty$のような取扱いに困る値が得られます．
物理的観察に基づく対策
小手先のテクニックで回避するのではなく，なぜこういうことが起こるのかを，物理的な観察，つまりその状況が表しているのはどういう現象か？を考えて，根本的な対策をする必要があります．
界面再構成関数と空間積分を再掲します．
f(x) = \frac{1}{2}\left\{ 1+ \tanh \left(\beta\left[X-d_i\right]\right)\right\}
\Phi_i=\int_{0}^{1}f(X)dX
Φ=1のとき
$\Phi_i=1$となるのは，$f(X)$が$X\in[0,1]$で常に$1$（つまり$\tanh(\beta[X-d_i])=1$）のときのみです．この状況を図に描いてみると，$d_i=-\infty$の意味が直ちに理解できます．
この図から，$d_i=-\infty$でなくても，浮動小数点数の有効桁数の範囲で$f(X)=1$が得られるような$d_i$の値で代用すればよいということが理解できます．
計算をしてみると，$\tanh(7\pi)$あたりで所望の値が得られることが分かります．余裕を見て，切りのよい数字$10\pi$を採用することにしましょう．
| y | tanh(y) |
|---|---|
| $ 1\pi$ | $0.99627207622074998$ |
| $ 2\pi$ | $0.99999302533961065$ |
| $ 3\pi$ | $0.99999998697517578$ |
| $ 4\pi$ | $0.99999999997567690$ |
| $ 5\pi$ | $0.99999999999995459$ |
| $ 6\pi$ | $0.99999999999999989$ |
| $ 7\pi$ | $ 1.0000000000000000$ |
| $ 8\pi$ | $ 1.0000000000000000$ |
| $ 9\pi$ | $ 1.0000000000000000$ |
| $10\pi$ | $ 1.0000000000000000$ |
界面再構成関数$f$の中では$y=(\beta\left[X-d_i\right])$ですが，$X$は$0$から$1$の範囲でしか変化しないので，影響が小さいとして無視すれば，
d_i = -\frac{10\pi}{\beta}
となります．つまり，$e_3-e_1 \approx 0$のときは，$\log$の計算を飛ばして$d_i=-\frac{10\pi}{\beta}$とすれば，$-\infty$を取り扱わなくてよくなります．
Φ=0のとき
$\Phi_i=0$となるのは，$f(X)$が$X\in[0,1]$で常に$0$（つまり$\tanh(\beta[X-d_i])=-1$）のときのみです．$\Phi_i=1$のとき同様にこの状況を図に描いてみると，$d_i=\infty$の意味が理解できます．同様に，$d_i$の値を$\infty$ではなく，浮動小数点数の有効桁数の範囲で$f(X)=0$が得られるような$d_i$の値で代用します．
$\tanh(-7\pi)$あたりで所望の値が得られるので，切りのよい数字$-10\pi$を採用します．
| y | tanh(y) |
|---|---|
| $- 1\pi$ | $-0.99627207622074998$ |
| $- 2\pi$ | $-0.99999302533961065$ |
| $- 3\pi$ | $-0.99999998697517578$ |
| $- 4\pi$ | $-0.99999999997567690$ |
| $- 5\pi$ | $-0.99999999999995459$ |
| $- 6\pi$ | $-0.99999999999999989$ |
| $- 7\pi$ | $- 1.0000000000000000$ |
| $- 8\pi$ | $- 1.0000000000000000$ |
| $- 9\pi$ | $- 1.0000000000000000$ |
| $-10\pi$ | $- 1.0000000000000000$ |
$X$の影響が小さいとして無視すれば，
d_i = \frac{10\pi}{\beta}
が得られます．$e_3e_1 \approx 1$のときは，$d_i=\frac{10\pi}{\beta}$とすれば，0除算を回避でき，かつ$\infty$の取扱いもなくなります．
まとめ
THINC法を用いるときは，界面位置$d_i$の値を$\pm\frac{7\pi}{\beta}$から$\pm\frac{10\pi}{\beta}$程度の範囲に抑えるようにするとよいでしょう．
よく分からない式でも，その式を使って表現したい物理的な現象は何かに着目すれば，解決策が見つかる場合もあります．