非ニュートン流体とは何ですか？

非ニュートン流体は、せん断応力とせん断速度の関係が非線形な流体です。見かけの粘度がせん断速度によって変化し、代表的な例として、かき混ぜると固くなるダイラタント性流体や、逆に流れやすくなる擬塑性流体（ケチャップなど）があります。CAE解析では、この複雑な粘性挙動を適切な構成則でモデル化する必要があります。

ニュートン流体と非ニュートン流体の違いは？

根本的な違いは、粘度が一定か変化するかです。ニュートン流体（水など）では粘度は定数で、せん断応力はせん断速度に比例します。非ニュートン流体では、見かけの粘度がせん断速度の関数となり、応力とひずみ速度の関係が非線形になります。この違いが流動挙動に大きな影響を与えます。

擬塑性流体（シアシニング）とは？

擬塑性流体は、せん断速度が増加すると見かけの粘度が低下する非ニュートン流体です。ケチャップや高分子溶液が例で、せん断が強まると高分子鎖が流れ方向に配向し、抵抗が減ることで起こります。この「せん断稀釈」現象は、CAEではパワーロー等のモデルを用いて数値シミュレーションされます。

非ニュートン流体のCAE解析で重要なことは？

連続の式や運動量保存則の形は同じですが、応力テンソルを定義する「構成則」が決定的に重要です。流体の種類（擬塑性、ダイラタント性など）に応じて、粘度がせん断速度にどう依存するかを適切な数式モデル（例：パワーローモデル）で表現し、運動方程式に組み込む必要があります。

非ニュートン流体 — CAE用語解説

カテゴリ: 用語集 | 2026-01-15

CAE visualization for non newtonian fluid - technical simulation diagram

非ニュートン流体

🧑‍🎓

先生、非ニュートン流体ってCFDでどう扱うんですか？

理論と物理

基本概念と支配方程式

🧑‍🎓

非ニュートン流体って、ニュートン流体と何が根本的に違うんですか？粘度が一定じゃないってことだけですか？

🎓

本質的な違いは、せん断応力とせん断速度の関係が線形か非線形かです。ニュートン流体では、粘度μは定数で、関係式は

$$ \tau = \mu \dot{\gamma} $$

です。一方、非ニュートン流体では、見かけの粘度ηはせん断速度

$$ \dot{\gamma} $$

の関数になり、

$$ \tau = \eta(\dot{\gamma}) \dot{\gamma} $$

と表されます。例えば、ケチャップは力を加えないと動かない（降伏応力あり）し、プールの水より水泳プールの方が抵抗が大きいのは、増粘剤で擬塑性流体にしているからです。

🧑‍🎓

「見かけの粘度」が変わるメカニズムは何ですか？分子レベルでは何が起きているんですか？

🎓

高分子溶液を例にとると、せん断速度が低い時は高分子鎖が絡み合って流れを大きく妨げます。せん断速度が上がると、鎖が流れ方向に引き伸ばされ、配向することで抵抗が減ります。これが擬塑性（シアシニング）です。逆に、懸濁液では粒子が高せん断領域で集まって「ジャミング」を起こし、抵抗が増すことがあります。これがダイラタント性です。数値的には、カラムーアモデルでは鎖の伸長を確率論的に扱います。

🧑‍🎓

支配方程式はナビエ-ストークス方程式そのままですか？それとも修正が必要なんですか？

🎓

連続の式と運動量保存則の形は同じですが、応力テンソルτの構成則（Constitutive Equation）が全く異なります。非圧縮性を仮定すると、運動方程式は

$$ \rho \left( \frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u} \right) = -\nabla p + \nabla \cdot \boldsymbol{\tau} + \rho \mathbf{g} $$

です。ここで、ニュートン流体では

$$ \boldsymbol{\tau} = 2 \mu \mathbf{D} $$

（Dは変形速度テンソル）ですが、非ニュートンではこの単純な関係が成り立ちません。この構成則をどうモデル化するかが全ての出発点です。

🧑‍🎓

構成則にはどんな有名なモデルがあるんですか？一番シンプルなのは？

🎓

パワーロー（べき乗則）モデルが最も基本的です。

$$ \eta(\dot{\gamma}) = K \dot{\gamma}^{n-1} $$

です。ここでKは粘度係数[Pa・s^n]、nはべき指数です。n<1なら擬塑性（シアシニング）、n>1ならダイラタント流体です。n=1ならニュートン流体に戻ります。例えば、0.3%のカルボキシメチルセルロース水溶液ではn=0.6程度、コーンスターチの濃厚懸濁液ではn>1のダイラタント性を示します。ただし、このモデルはせん断速度が極端に低いor高い領域では破綻します。

数値解法と実装

離散化とソルバー設定

🧑‍🎓

非ニュートン流体のCFD解析では、粘度が場所によって変わるので、離散化で特別な処理が必要ですか？

🎓

はい、重要な点は「見かけの粘度」が速度勾配（せん断速度）の関数であるため、非線形性が強くなることです。有限体積法では、セル界面での粘度η_faceをどう補間するかが鍵です。単純な算術平均では不正確で、せん断速度も界面で再計算するか、アップウィンド的な考え方を取り入れることがあります。また、せん断速度

$$ \dot{\gamma} = \sqrt{2 \mathbf{D}:\mathbf{D}} $$

の計算には、速度勾配テンソルの不変量を用いますが、これもセル中心と界面で一貫性を持たせる必要があります。

🧑‍🎓

ソルバーの収束が悪くなりそうですが、ニュートン法以外に有効な解法はありますか？

🎓

擬似時間進行法（Pseudo-transient）が有効です。定常解析でも、非物理的なダンピング項を加えた時間発展方程式を解き、定常状態に収束させます。OpenFOAMの`simpleFoam`非ニュートン対応ソルバーでは、粘度の更新を反復の外側ループ（outer iteration）で行い、内側の速度-圧力の反復（SIMPLE）では粘度を固定する「分割解法」を採用しています。また、べき乗則モデルでnが0.3など非常に小さい場合、初期条件をニュートン流体（n=1）からゆっくりnを下げていく「継続法」を使うこともあります。

🧑‍🎓

降伏応力があるビンガム流体みたいなモデルでは、未降伏領域（固体的挙動）と降伏領域（流動）の境界をどう扱うんですか？

🎓

それが一番難しい点です。ビンガムモデルでは、応力が降伏応力τ_yを超えた領域だけが流動します。数値的には、非常に高い粘度を持つ疑似固体領域として扱う「粘度正則化」が一般的です。例えば、Papanastasiouの正則化モデル：

$$ \eta(\dot{\gamma}) = \mu_p + \frac{\tau_y}{\dot{\gamma}} \left[ 1 - \exp(-m \dot{\gamma}) \right] $$

を使います。パラメータmが大きいほど理想的なビンガム流体に近づきますが、数値的剛性も増します。m=1000 1/s程度から始めるのが現実的です。

実践ガイド

ワークフローとチェックリスト

🧑‍🎓

実際に非ニュートン流体の解析を始めるとき、最初に何を決めるべきですか？

🎓

絶対に必要なのは「実測の流動曲線データ」です。せん断速度の範囲を明確にしましょう。例えば、円管内流れの壁面せん断速度は

$$ \dot{\gamma}_w = \frac{8V}{D} $$

（Vは平均流速）で概算できます。レオメーターで測定したτ-γ˙データを、パワーロー、キャリオー、ヘルシェル・バルクレーなどのモデルにフィッティングします。フィッティングの決定係数R²が0.99を切るようなら、モデルの適用範囲を疑うか、複数のモデルを領域ごとに切り替える必要があります。

🧑‍🎓

メッシュの切り方に注意点はありますか？特に壁面近くは？

🎓

壁面せん断速度が最も高くなるため、粘度が大きく変化する境界層を捉える必要があります。ニュートン流体のy+<1のような指針はそのまま使えません。まずはニュートン流体としてシミュレーションし、得られたせん断速度分布を見て、粘度が急変する領域（例えば、ビンガム流体なら降伏面付近）にメッシュを集中させます。特に、円管中心でせん断速度が0になる場合、そこで粘度が発散しないよう、モデルの正則化が機能するか確認する必要があります。

🧑‍🎓

結果の検証はどうすれば信頼できますか？実験データがなければ？

🎓

実験がなくても、理論解や既知の相関式と比較する「ベンチマーク」は必須です。例えば、円管内のパワーロー流体の発達流れ圧力損失は、理論的に計算できます（レイノルズ数は一般化レイノルズ数Re_gを使用）。また、メッシュ依存性調査では、圧力損失や中心流速など重要な出力量がメッシュ数を2倍にしても変化が1%以内になることを確認します。ソルバーの残差が下がっても、グローバルな力のバランス（流入運動量と圧力勾配+壁面せん断力）がとれているかチェックしてください。

ソフトウェア比較

Ansys, Abaqus, COMSOL等の実装

🧑‍🎓

Ansys Fluentで非ニュートン流体を扱う場合、どのモデルが標準で使えますか？

🎓

Fluentには豊富な非ニュートンモデルが組み込まれています。主なものは、パワーロー、キャリオー、ヘルシェル・バルクレー（降伏応力付き）、クロス、バード・キャリオーなどです。ユーザー定義関数（UDF）で独自モデルも実装可能です。実務では、高分子溶融樹脂の射出成形シミュレーションでは「Modified Cross Model」がよく使われます。設定時は、せん断速度の下限・上限にクリップ値を設定する（例えば1e-6 1/sと1e6 1/s）ことで、数値的发散を防ぎます。ソルバー設定は「Pressure-Based」が一般的で、勾配の再構成スキームは「Least Squares Cell Based」より「Green-Gauss Node Based」の方が安定しやすいです。

🧑‍🎓

COMSOL Multiphysicsは「マルチフィジックス」が売りですが、非ニュートン流体連成では何が得意なんですか？

🎓

COMSOLの強みは、非等温による粘度変化（温度依存性）や化学反応による組成変化を、流体流れと完全に連成して解ける点です。「非等温非ニュートン流れ」インターフェースでは、粘度モデルとしてArrhenius型の温度依存性

$$ \eta(T) = A \exp(E_a / RT) $$

を、せん断速度依存性と組み合わせて定義できます。また、自由表面を伴う塗布プロセスでは、レベルセット法やフェーズフィールド法と連成させ、表面張力と非ニュートン性を同時に考慮した解析が可能です。ユーザー定義の粘度関数をMATLABライブリンクで直接記述できるのも特徴です。

🧑‍🎓

構造解析のAbaqusで非ニュートン流体を扱うことはあるんですか？

🎓

はい、特に「キャビティ充填解析」や「ゴム材料の超弾性+粘弾性解析」で重要です。Abaqus/CFDモジュールもありますが、StandardやExplicitの「粘性」材料モデルとして定義する場合が多いです。例えば、樹脂の射出成形シミュレーション（Abaqus/ExplicitのCEL法を用いる場合）では、7係数のCross-WLFモデルが業界標準です。このモデルはせん断速度と温度の両方に依存します：

$$ \eta(\dot{\gamma}, T) = \frac{\eta_0(T)}{1+(\eta_0 \dot{\gamma} / \tau^*)^{1-n}} $$

ここでη_0はゼロせん断粘度で、WLF方程式で温度依存性を表現します。材料データはMoldflowなどの専用ソフトからインポートすることもあります。

トラブルシューティング

よくあるエラーと対策

🧑‍🎓

解析中に「負のせん断速度が検出されました」というエラーが出ます。せん断速度は本来負にならないはずでは？

🎓

数値計算上の丸め誤差や、発散途中の非物理的な速度場によって、セル内の変形速度テンソルDの第2不変量の計算結果が微小な負の値になることがあります。対策は二つあります。1) せん断速度の計算式に小さな正の数εを加えてクリップする：

$$ \dot{\gamma} = \sqrt{2\mathbf{D}:\mathbf{D} + \epsilon^2} $$

(ε=1e-6 1/sなど)。2) より根本的には、初期条件や境界条件を見直し、流速が急激に変化する場所のメッシュを細かくして、速度勾配の計算精度を上げます。Fluentでは「Shear Rate Limit」を設定するのが簡単です。

🧑‍🎓

べき乗則モデルでnが小さい（0.2など）と、ソルバーの残差が振動して全く収束しません。どうすればいいですか？

🎓

nが小さいと、せん断速度が低い領域で粘度が非常に高くなり、方程式の剛性が極端に高まることが原因です。対策ステップ：1) **初期化**：まずn=1（ニュートン流体）で完全に収束する解を得る。2) **継続法**：解を初期値として、nを0.9, 0.8,...と段階的に減らしながら解析を連続して実行する。多くのソルバーで「パラメトリックスタディ」機能を使えば自動化できます。3) **緩和係数**：速度と圧力のアンダーリラクゼーション係数を大幅に下げる（0.2以下）。4) **モデル変更**：現実的には、低せん断域で粘度が有限値に収束する「キャリオーモデル」などの方が適している可能性があります。

🧑‍🎓

ビンガム流体の解析で、降伏していないはずの領域（プラグ領域）が数値的にわずかに流れています。これは許容誤差ですか？

🎓

それは「正則化誤差」です。先述のPapanastasiouモデルでは、パラメータmが有限であるため、せん断速度がゼロに近い領域でも粘度が無限大にならず、ごく微小な流れ（疑似クリーップ）が生じます。この誤差を評価するには、プラグ領域の最大流速V_plugを代表流速Vで無次元化し、

$$ V_plug / V < 0.01 $$

程度か確認します。許容できない場合は、正則化パラメータmを大きくしますが、その分ソルバーの収束は悪化します。実用的には、この微小流れが工学的結論（圧力損失、混合効率など）に影響を与えないかどうかで判断します。

🧑‍🎓

ユーザー定義粘度関数を実装したら、計算コストが異常に高くなりました。原因は？

🎓

主な原因は二つ考えられます。1) **関数内での条件分岐（if文）**：セルごとに「せん断速度が閾値以上か」などで分岐すると、ベクトル化が阻害され、CPUキャッシュの効率が大幅に低下します。可能な限り、滑らかな関数（tanhなど）で近似して分岐をなくします。2) **高次の数学関数の多用**：exp(), log(), pow() は計算コストが高いです。特にpow(x, n)でnが非整数の場合、内部でexp(log(x)*n)が実行されます。事前にテーブル参照（ルックアップテーブル）を用意し、補間で粘度を求める方法に切り替えると、速度が10倍以上向上することもあります。Ansys FluentのUDFなら`DEFINE_PROPERTY`マクロで、`real shear_rate = sqrt(2.0*C_DOT(c,t)*C_DOT(c,t));`のようにせん断速度を取得しますが、この計算自体もコストです。