非ニュートン流体 — 数値解法と実装

根本的な違いは、粘度が速度場に依存することだ。Navier-Stokes方程式の粘性項が非線形になる。

$\mathbf{D} = \frac{1}{2}(\nabla\mathbf{u} + \nabla\mathbf{u}^T)$ はひずみ速度テンソルだ。

解法としてはPicard反復（逐次代入法）が基本になる。

1. 前のステップの速度場から $\dot{\gamma}$ を計算

2. $\eta(\dot{\gamma})$ を更新

3. 更新された粘度でN-S方程式を解いて新しい速度場を得る

4. 収束するまで1-3を繰り返す

そうだ。この外側の反復ループが収束しにくいのが非ニュートン流体計算の難しさだ。特にshear-thinning度が強い（$n$ が小さい）場合や降伏応力がある場合は要注意。

非ニュートン流体ではレイノルズ数の定義も変わる。Power-law流体の管内流ではMetzner-Reed一般化レイノルズ数が使われる。

層流-乱流遷移は $\text{Re}_{\text{MR}} \approx 2100$ で起きる（ニュートン流体とほぼ同じ臨界値）。ただし遷移後の乱流挙動はニュートン流体とは大きく異なる。

CFDで非ニュートン流体を解く際の離散化の注意点を説明しよう。

粘度のフェース値補間が重要になる。セル中心で計算した粘度をフェースに補間する際:

• 調和平均: 粘度が急変する界面で推奨。$\eta_f = \frac{2\eta_L \eta_R}{\eta_L + \eta_R}$
• 算術平均: 粘度変化が緩やかな場合。$\eta_f = \frac{\eta_L + \eta_R}{2}$

shear-thinning流体では、壁面近傍と流路中心で粘度が数桁異なることがある。例えば高分子メルトの射出成形では $\eta$ が $10^1 \sim 10^4\,\text{Pa}\cdot\text{s}$ の範囲で変化する。この急変に対応するため、壁面近傍のメッシュを十分に細かくする必要がある。

その通り。緩和係数を0.3-0.5に下げ、粘度の更新にも緩和を入れるのが実務的なテクニックだ。