層流

目次	粘弾性流体の流れ
はじめに	エネルギー輸送 - 非等温流れ
概要的な例と解析タイプ	流れ境界条件
方程式	流体流れモデルの収束
ニュートン流れと非ニュートン流れ	参考文献

はじめに

流体の解析と挙動は科学と工学において基本的な重要性がある．このモノグラフでは偏微分方程式を使った流体のモデル化の基本を紹介する．流体力学解析を行うのに必要な方程式および境界条件を導出し説明する．

流体はそれぞれ空気や水のような気相および液相の物質である．「流体」という言葉は液体と気体の両方を包含する．流体の反対は固体であり，その違いを示す特性は流体は静止時にせん断力に抵抗できないということである．流体は容器の形になる．流体と固体の間の境界線は必ずしもはっきりとはしていない．例えば，アスファルトは通常固体と考えられるが，長時間経つと流体のように振舞い始める．

通常以下のように流体流れを分ける：

クリープ流

層流

乱流

これら3つの流れの様式の違いは，慣性力と粘性力の比であるレイノルズ数で与えられる．これは流体流れに乱流が生じるかどうかを予測するために使われる．乱流は大きいレイノルズ数で特徴づけられ，層流はレイノルズ数が小さいときに発生する．クリープ流という言葉は非常に小さいレイノルズ数を持つ流れに使われる．レイノルズ数の概念は後でより詳しく説明する．また，レイノルズ数が大きいまたは小さいということは何を意味するのかについても触れる．

このモノグラフでは層流に焦点を当てる．

偏微分方程式(PDE)を使って流体をモデル化することだけが流体のモデル化の方法ではない．他の方法として常微分方程式(ODE)を設定するというものがある．Wolfram System Modelerではこのアプローチに従っている．大雑把に言うとSystem Modelerのアプローチは管の中の流れのような大きな流れの系の相互作用により適しており，偏微分方程式のアプローチは特定の流体流れ挙動の細かい解析に適している．両方のアプローチを組み合せて使った方がよい場合もある．

ここで取るアプローチは，入門セクションでは単独の流体流れであるドリブンキャビティ流れの例を使って，利用できるさまざまな流体流れの解析タイプ機能を紹介する．その後で根底にあるアイディアと概念をより理論的に説明する．理論的背景は，さまざまな解析タイプに対する直感ができるとより理解しやすくなる．それから利用可能な境界条件を説明する．

流体解析の目標は，制約条件下の流体のベクトル値の流速を求めることである．他に，非圧縮性流体と圧縮性流体のどちらを考えるかによって，それぞれスカラー圧力場やスカラー密度場を求めるということもある．流体挙動の解析と解釈は問題となる流れ挙動についてよりよい品質の工学設計を行うのに役立つ．例えば物体のある部分では十分な流れがないことがある．このような部分を特定し改善することができる．

流体解析は通常段階を追って行われる．まず，流体が流れる物体では，幾何学的モデルを作成する必要がある．幾何学的モデルは通常コンピュータ支援設計(CAD)プロセスで作成される．CADモデルはインポートすることも製品内で作成することもできる．形状をインポートする場合，DXF，STL，STEP等の一般的なファイル形式がサポートされている．これらの形状はImportでインポートすることができる．他の方法として，例えばOpenCascadeLinkを使って製品内で幾何学的モデルを作成するというものがある．幾何学的モデルが使えるようになると，どのような解析を行うかを考える必要がある．現在サポートされている解析タイプは静的で時間依存の流体解析である．次のステップは境界条件と制約条件の設定である．問題となっている流体の材料特性によってさらにPDEモデルが指定される．PDEモデルが完全に指定されると，それに続く有限要素法が調査中の所望の数量を計算する．その後この数量は後処理される．可視化されることもあれば導出される数量が計算されることもある．このノートブックはCADモデル生成以外のすべての必要なステップを示す．

そのようなモデリングプロセスは最終的にNDSolveおよびParametricNDSolveで解くことのできるPDE系になる．

このノートブックで示すシミュレーション結果のアニメーションの多くはRasterizeを呼び出して生成されている．これによってこのノートブックが必要とするディスク容量を減らしている．欠点は，アニメーションの視覚的品質が呼び出さない場合よりも劣るということである．高品質のグラフィックスが入手したい場合は，Rasterizeの呼出しを削除するかコメントアウトするかするとよい．

高忠実度の可視化を得るためには，ラスタライズの処理をコメントアウトする．

このチュートリアル全体で使われている記号とそれに対応する単位は用語集セクションにまとめてある．

有限要素パッケージをロードして$HistoryLengthを設定する：

高忠実度の可視化を得るためには，ラスタライズの処理をコメントアウトする．

このチュートリアル全体で使われている記号とそれに対応する単位は用語集セクションにまとめてある．

概要的な例と解析タイプ

このモノグラフは，実際の方程式や理論よりもワークフローを重視する概要的な例から始める．この例では細い管の中のグリセロールの流れを考える．

流体モデルの作成には常に同じステップが必要である：

形状の作成

材料の選択

境界条件の指定

メッシュ生成

偏微分方程式を解く

解の後処理

形状を設定して可視化する：

モデル変数を設定する：

従属速度変数，はそれぞれ方向，方向の変位を表す．は圧力である．

流体はグリセロールである：

ここで化学実体を使う．その代りに材料特性を指定することもできる．

材料特性から流体パラメータを指定する：

流体流れPDE成分を設定する：

方向の流入速度は境界条件で規定される．これは従属変数の流入プロファイルを領域の左側に指定することで行える．領域の下と上の部分には壁面境界条件を指定する．これを行うためには，これらの部分のおよび成分を0に設定する．領域の右側には流出条件を設定する．ここで圧力は任意に0に設定する．一般にそれぞれの従属変数にはディリクレ条件かロビン型ノイマンを指定する．ノイマン0のような純粋なノイマン値では不十分である．この場合方程式の解は定数までしか正しくない．

境界条件を設定する：

速度が圧力よりも高次で補間されると安定解が見付かる．NDSolveを使うと，それぞれの従属変数の補間次数が指定できる．ほとんどの場合"InterpolationOrder"の指定は必要ではないが，流体流れに有限要素法を使う場合はこれが標準的な方法であり，Wolfram言語ではこのように行う．この設定はP2P1型またはQ2Q1型の要素と呼ばれるものに対応しており，これらの要素はTaylor–Hood（テイラー・フード）要素としても知られる．"InterpolationOrder"を使わないと，すべての変数に対して二次補間を使うことになり，不安定性につながる．

PDEを解く：

多くの場合メッシュを生成する必要はなく，形状をソルバに渡すだけで十分である．

流体速度を可視化する：

圧力分布の等高線を可視化する：

これがワークフローの概要である．次のセクションではそれぞれのステップをもう少し詳しく説明する．

定常解析

このドリブンキャビティの例では，流体力学の古典的なベンチマークを調べる[6]．矩形領域が流体で満たされている．上部にはコンベアーベルトのように，正の方向に一定の流速で流体を流す装置がある．残りの辺は壁面であり滑りなしの条件を持つ．滑りなしとは壁面協会の速度が方向と方向の両方で0ということである．左下角には基準圧条件がある．

流体で満たされた矩形領域で，上部は方向に速度で流れる．他のすべての辺は壁面であり，滑りなし境界条件()が設定されている．

このドリブンキャビティは時間非依存の定常モデルとしてモデル化される．従属速度変数と，従属圧力変数を設定する．独立変数は空間変数のとである．

定常モデル変数を設定する：

矩形領域を作成する：

標準的なテキストではこの例をさまざまなレイノルズ数を使って提示している．他のすべてが同じままだと，レイノルズ数の値だけによって流れのパターンが変化する．

の流体流れパラメータ pars を設定する：

流体流れPDEを設定する：

上部の境界条件を設定する：

壁面境界条件を設定する：

基準圧点は「穏やかな」部分，つまりおもしろい挙動から離れた部分に適用する．

圧力点境界条件：

境界条件を組み合せる：

流体流れPDEを解く：

返された解は，方向の速度プロファイルである速度プロファイル，方向の速度プロファイルである速度プロファイル，圧力分布を制約する．解を可視化する方法はいくつかある．

NDSolveに与える"InterpolationOrder"メソッドオプションは，両方の速度が2次で補間され圧力が1次で補間されることを指定する．"InterpolationOrder"オプションの指定は，有限要素でテイラー・フード要素として知られるものと同等であり，解の安定化に役立つ．残念ながら特定のPDEに対してこのオプションを自動的に設定する方法はないので，FluidFlowPDEComponentのようにナビエストークス方程式に基づく流体流れPDEについてはこのオプションを設定するようお勧めする．

後処理

一般的な可視化テクニック

計算された圧力はスカラー場なので，他のスカラー場のように可視化することができる．これを可視化するには等高線プロットまたは密度プロットが向いている．

等高線プロットで圧力場を可視化する：

速度場はベクトル値場である．ここではベクトルプロット，流線プロット，線積分畳み込みプロットが可視化に向いている．

速度場をベクトルプロットとして可視化する：

速度場を流線プロットとして可視化する：

速度場を線積分畳み込みプロットとして可視化する：

二次元では速度プロットは便利な情報も提供する．

流れ場の渦度を計算する関数を作成する：

渦度の[-10,10]の範囲で等高線を可視化する：

パラメトリック解析

次のステップとして，さまざまなレイノルズ数に対する流れ場を比較する．ParametricNDSolveはNDSolveと同じように使えるが，パラメータを指定することができる．パラメータに数値を入れると，ParametricNDSolveは数値解を計算する関数を返す．この場合はレイノルズ数をパラメータとして使う．

記号的レイノルズ数 re を設定する：

流体流れPDEを設定する：

さまざまなレイノルズ数に対するParametricFunctionを設定する：

いくつかのレイノルズ数の解を計算する：

NDSolveを使った前の例では，{uVel,vVel,pressure}の形式のリストに解を保存した．ここではすべてのシミュレーションの解を変数 solutions に保存する．解のベクトルのそれぞれの項目には場，場，圧力の補間関数が含まれている．

それぞれの解の速度場を可視化する：

レイノルズ数が大きすぎると，PDEが非常に強非線形になるため，非線形PDEの収束が妨げられる．FindRootはNDSolve内で非線形方程式を最終的に解く関数なので，収束の失敗は通常FindRootからのメッセージで知らされる．しかし望みがないわけではない．モデルが収束する範囲を拡張するために取る手段がいくつかある．

すべてのレイノルズ数がすぐに収束する訳ではない：

収束の失敗を回避する方法はいろいろあり，セクション「流体流れモデルの収束」で説明してある．ここでは解析に使われるメッシュを詳しく見てみる．上の例ではメッシュは自動的に生成された．一つの方法は，手元の問題により適したメッシュを使うことである．「より適した」というのは手元の問題によって異なる．この例では速度場と速度場の勾配は上部と右側壁面で最も大きい．

最後に計算されたレイノルズ数に対する

速度場と

速度場の最大勾配の位置を求める：

ここで生成するメッシュは，急勾配においてより細かい解像度を持つことによって急勾配を反映する．流れ勾配が壁面で高くなるのを見ると，それを向上する可能な方法はグレーデッドメッシュを使うことである．これは関数ToGradedMeshを使って行うことができる．これにより境界での微調整をしながら全体的な要素数をそれほど増やさないでおける．MeshRefinementFunctionを使う等の別の調整方法も使える．メッシュ生成処理についての詳細は「要素メッシュの生成」チュートリアルに記載されている．

壁面と上部を調整したメッシュを生成する：

メッシュを可視化する：

勾配メッシュでは，全体的な要素数を減らすと同時に壁面付近の要素密度を増加させた．

構築されたメッシュを使って新しいパラメータ関数を作成する：

さまざまなレイノルズ数に対する解を計算する：

StreamPointsとして使われるいくつかのメッシュ点を抽出する：

さまざまなレイノルズ数における解を可視化する：

アニメーションの品質を向上させる方法についてはこちらをご覧いただきたい．

ドリブンキャビティは一般的な流体力学のベンチマーク問題[6]なので，数値シミュレーションを比較する文献データが存在する．

ベンチマーク結果をインポートする：

において計算された解とさまざまレイノルズ数についてのベンチマーク結果を比較する：

において計算された解ととさまざまレイノルズ数についてのベンチマーク結果を比較する：

解は文献の値と一致する．

レイノルズ数10000についての右下の速度場を調べる：

レイノルズ数10000についての左下の速度場を調べる：

レイノルズ数10000についての左上の速度場を調べる：

ドリブンキャビティの例におけるより大きいレイノルズ数の計算は「流体流れモデルの収束」セクションを参照のこと．

時間依存解析

次のワークフローの例として，時間依存のドリブンキャビティの例を解く．ここでは流体は最初静止しており，駆動装置は時間とともに強化される．その他はすべて同じままである．

時間依存変数を次のように設定する：

過渡ナビエストークス演算子をレイノルズ数1000について設定する：

流入境界条件が初期条件と一貫するようにするために，時間経過とともに流入をスムーズに増加させるヘルパー関数を作成する．

時間とともに流入境界条件を増大させる関数：

上部の境界条件を設定する：

壁面境界条件を設定する：

圧力点境界条件：

境界条件を組み合せる：

最初流体は静止している．初期条件は0に設定する：

時間積分を監視することができる．

ナビエストークス方程式を時間積分する：

さまざまな時間についての速度場を可視化する：

アニメーションの品質を向上させる方法についてはこちらを参照されたい．

方程式

流体流れのシミュレーションを行う場合の最初のステップは，どのフローレジームが想定されるかを分類することである．このために流れ挙動のさまざまな面を捉える特性数を使う．一つの特性はレイノルズ数で行われる．レイノルズ数は流体の内部力の粘度に対する割合である．

ここで []は動粘度， []は密度， []は典型的な流速， []は代表長さである．レイノルズ数はおおよそどのようなタイプの流れが想定できるかを示してくれる．

の流れはクリープ流れレジームである．レイノルズ数が等と大きい場合の流れでは，乱流となり得る．どのレイノルズ数で流れが層流から乱流に変化するかは，その流れの指定による．例えば管の中ではで乱流になる可能性がある．それとは逆に平板上の流れはで乱流になることがある．

通常流体のフローレジームは以下のように分ける：

クリープ流：

層流：

乱流：

水と比べ，蜂蜜は粘性流体である．流体の粘性が大きいほど，層流挙動を示す可能性が高い．一方気体は粘性が低いためレイノルズ数は無限に近付く．したがって，ほとんどの分野では気体は非粘性と考えられている．つまり簡単に言うと気体の粘性は0と想定される．気体は非粘性流体とも呼ばれる．しかしこのような想定が不可能な分野もある．

粘性層流はナビエストークス方程式というベクトル値の方程式の組で記述される．広範な流体流れ現象は，層流以上でさえもナビエストークス方程式で記述することができる．ナビエストークス方程式は流体力学において包括的な方程式である．ナビエストークス方程式は偏微分方程式の集合であり，関心の対象となっているフローレジームによってさまざまな形式で存在する．

非粘性流れはオイラー方程式で記述することができるが，ここでは扱わない．

ナビエストークス方程式は運動量方程式と連続方程式の組合せである．

圧縮性流体流れの運動方程式は粘性応力テンソル（単位は[]）で表されるベクトル値の方程式である：

ここで

は単位[]の密度

は単位[]の流体流れ速度ベクトル

は単位[]の時間

は単位[]の粘性応力テンソル

は単位[]の流体圧力

は単位[]の体積力密度

は単位行列

である．粘性応力テンソルは以下で与えられる：

ここで

は単位[]の動粘度

は単位[]のひずみ速度テンソル

である．は以下で与えられる：

粘性応力テンソルの定義を運動方程式に挿入すると以下が得られる：

ナビエストークス運動方程式の形式の中には以下のような第二粘度または体積粘度と呼ばれる追加の項を含む場合がある：

しかし，ほとんどの実際の目的ではは無視でき，である．

運動方程式の単位は，単位[]の力密度である．運動保存は実質的に流体についての運動の第二法則である．これは方程式を少し並べ替えると分かる：

運動方程式はほとんどの場合，質量保存を記述する質量連続方程式と一緒に解く：

ここで密度は従属変数であり流体流れ速度ベクトルはである．

ナビエストークス方程式は，領域の最小の幾何学的特徴が流体の分子の平均自由行程よりずっと大きい場合のみ有効である．

ナビエストークス方程式は圧縮性，弱圧縮性，非圧縮性の流れに使える．弱圧縮性の流れとは，密度に対する圧力の依存性が無視でき，基準圧力で評価されるものである．

質量密度の定数値では，質量連続方程式は体積連則方程式に簡約され，それによって粘性応力テンソルはに簡約される．が定数のとき，発散はなのでをと書くことができる．

ナビエストークス方程式の非圧縮性形式は以下である：

成分形式で書いた非圧縮性ナビエストークス方程式は以下で与えられる：

一定の動粘度のとき，さらに以下のように簡約できる：

ニュートン流れと非ニュートン流れ

非ニュートン流体の粘性は適用される力で変化する．水や蜂蜜等のニュートン流体の場合，粘性は水をどのくらいかき混ぜるかによって粘性は変化しない．非ニュートン流体では，攪拌の量で粘性が変化する．材料によっては混ぜられると混ぜにくくなるものもあれば混ぜやすくなるものもある．実際すべての気体およびほとんどの液体はニュートニアンと考えられる．非ニュートン流体にはペンキ，血液，水とでんぷんの混合液，液体ポリマー等がある．

次のセクションでは，さまざまな非ニュートン流体流れモデルを紹介する．まず圧縮性ナビエストークス方程式の変形を始める：

これで因子を導入する：

ここで単位[1/s]のひずみ速度テンソルを定義することができる：

このひずみ速度は固体力学で定義した微小ひずみ測定によく似ている．固体力学では従属変数は単位[m]の変位であり，は微小ひずみを定義する．流体力学では従属変数は単位[m/s]の速度，はひずみ速度テンソルである．ひずみ速度テンソルは速度テンソル，ひずみテンソルと言われることもある．

ニュートン流体の粘性応力テンソル（単位[Pa]）は以下である：

これにより異なる形式のナビエストークス方程式が与えられる：

非ニュートン流体では応力ひずみ速度関係は非線形である．この非線形性は見かけ粘度で表すことができる．見かけ粘度はせん断速度（単位は[1/s]）の関数である．見かけ粘度は実効粘度とも呼ばれる．

簡単にするために，一般性についての制約することなく，項のない非圧縮性非ニュートン流体を仮定すると以下が得られる：

せん断速度は次で計算される：

ここでは収縮であり以下で計算される：

せん断速度を計算する関数は次のように実装され，PDEModels`コンテキストの中にある：

ShearRate[strainRate_] := 
    Sqrt[2 * TensorContract[strainRate * strainRate, {{1}, {2}}]];

動粘度は温度と圧力に依存することがあり，非ニュートン流体では動粘度はさらにせん断速度にも依存する．

さまざまな非ニュートン構成モデルを使って見かけ粘度を表す．以下の非ニュートン流体流れモデルを実装する：

ベキ乗則

Carreau

Bingham-Papanastasiou

Herschel-Bulkley-Papanastasiou

カスタムの見かけ粘度関数

カスタムの粘性応力テンソル関数

以下で示すモデルにはすべてパラメータ pars に"FluidDynamicsMaterialModel"の項目が必要である．

ベキ乗則

ベキ乗則モデルはよく使われる非ニュートン流体モデルである．

ベキ乗則は次を実装する：

ここではデフォルトの動粘度になる因数である．基準せん断速度はデフォルトの1 [1/s]である．デフォルトの最小せん断速度は [1/s]に設定される．スカラーは指数である．のとき，ニュートン流体が得られベキ乗則が一般化ニュートン流体モデルを実装する．のときはせん断減粘流体（擬塑性流体とも呼ばれる）となる．ペンキ，血液，ケチャップはこのカテゴリに属する．のときは，せん断増粘流体（ダイラタント流体とも呼ばれる）となる．せん断増粘流体の例としては，水とコーンスターチの混合液や流砂がある．

ここで使うモデルは文献で通常みられる元のは少し異なっている．もとのモデルはせん断速度がゼロになるときの粘性無限大を予測する．実際の流体では，せん断速度は一定値になりやすい．もとの式のこの欠点を避けるための式がここで与えたものである．ベキ乗則の欠点はニュートンのプラトーが記述できないことである．Carreauモデルではできる．

ベキ乗則モデルのパラメータ名：

ベキ乗則は動粘性率の式として与えられることがよくある：

ここでは整合度と呼ばれる．

せん断流れとニュートン流れの図を以下に示す．

ニュートン流体のせん断増粘流れとせん断減粘流れについてのひずみ速度とせん断応力の比較．

せん断減粘流体では，ひずみ速度が増大すると見かけ粘度が減少し，せん断増粘流体ではひずみ速度の増大で見かけ粘度が増大することが分かる．

非ニュートン流体に対するベキ乗則モデルの実際の実装は以下に示す通りである．このモデルの構築方法は，他のニュートン流体モデルの構築の例として使うことができる．

ベキ乗則非ニュートン流体モデルの実装：

ApparentViscosity["PowerLaw", vars_, pars_, data__]:=                           
Module[                                                                         
    {model, strainRate, shearRate, m, apparentViscosity, minShearRate, n, refShearRate},
    model = pars["FluidDynamicsMaterialModel"]["PowerLaw"];
    minShearRate = model["MinimalShearRate"];
    refShearRate = model["ReferenceShearRate"];
    n = model["Exponent"];
    m = model["PowerLawViscosity"];
    strainRate = "StrainRate" /. data;
    shearRate = PDEModels`ShearRate[strainRate];
    apparentViscosity = m * (Max[shearRate, minShearRate]/refShearRate)^(n-1);
    apparentViscosity
]

ベキ乗則の例

次のセクションではベキ乗則の非ニュートン流体の例を与える．

拡大領域を設定する：

領域をメッシュにし，可視化する：

ここの例題では動粘度が指定されている．密度が1として与えられているので，動粘度は粘度に等しい．

実験するために，ベキ乗則指数およびベキ乗則因数を記号パラメータにする．

変数，パラメータ，パラメトリックベキ乗則PDEを設定する：

流入プロファイルを{1/2,0}として与え，流出での圧力を0に設定する．壁面は滑りなしとする．

境界条件を設定する：

パラメトリック関数を作成する：

さまざまな指数

と因数

についてパラメトリックモデルを評価し，かかる時間を測定する：

x方向の速度を抽出する：

チャンネルの中央から上までの流束プロファイルを1にスケールしてプロットする：

Carreau

Carreauモデルも一般化ニュートン流れモデルを実装する．これはポリマー流れをモデル化するのに便利である．Carreauモデルは次を実装する：

ここでは無限せん断速度粘度，は零せん断速度粘度，は緩和時間（単位は[s]），はベキ乗則の指数，はニュートンモデルとベキ乗則モデルの間の移行を制御する移行指数である．またはのとき，このモデルはニュートン流体の標準モデルに還元される．ではCarreauモデルはニュートン流れに還元される．ではCarreauモデルはベキ乗則モデルに近付く．

ここで示すモデルは，通常Carreau-Yasudaモデルと呼ばれるCarreauモデルの一般化である．もとのCarreauモデルはデフォルト設定ので得られる．

Carreauモデルを使うと非常に小さいせん断速度についても非常に大きいせん断速度についてもニュートンプラトーの記述が可能である．ニュートンプラトーは非ニュートン流体と関連のある現象である．このような流体は非常に小さいせん断速度では一定粘性を示す．これは第一プラトーと言われる．第二プラトーはせん断速度が非常に大きい場合に見られる．

Carreauモデルのパラメータ名：

"TransitionExponent"パラメータはオプショナルである．これが与えられていない場合は，デフォルト値のが使われる．

非ニュートン流れに対する実際のCarreauモデルの実装は以下に示す．このモデルの構築方法は，他の非ニュートン流れモデルの構築の例として使うことができる．

Carreau非ニュートン流れモデルの実装：

ApparentViscosity["Carreau", vars_, pars_, data__] :=
 Module[
      {model, mu0, muInf, n, a, lambda, strainRate, shearRate, apparentViscosity},
      model = pars["FluidDynamicsMaterialModel"]["Carreau"];
      mu0 = model["ZeroShearRateViscosity"];
      muInf = model["InfiniteShearRateViscosity"];
      n = model["PowerLawExponent"];
      a = model["TransitionExponent"];
      lambda = model["Lambda"];
      strainRate = "StrainRate" /. data;
      shearRate = ShearRate[strainRate + $MachineEpsilon];
      apparentViscosity = muInf + (mu0 - muInf) *
             (1 + (lambda * shearRate)^a )^((n - 1)/a);
      apparentViscosity
  ]

モデルのパラメータを詳しく見ていこう．以下ではパラメータの選択がどのように見かけの粘度に影響を及ぼすかを示すいくつかのプロットを提示する．最初に，固定されたデフォルトのパラメータを持つモデルを考えてから，パラメータを1つずつ変化させ，その結果の見かけの粘度を比較する．これにより，パラメータの選択がどのようにモデルに影響を与えるかが分かる．パラメータのデフォルトとして以下の値を選ぶ：

見かけの粘度の関数を定義する：

最初に変化させるパラメータはである．このパラメータは，せん断速度がゼロのときに流体の粘度を記述するために使われる．せん断速度がゼロということは流体の速度が一定ということである．

基準値を設定して，パラメータ値

をテストする：

[0,3]の区間での

のさまざまな値に対する見かけの粘度を比較するプロットを作成する：

このプロットから，小さいせん断速度のときが見かけの粘度に影響を及ぼしていることが分かる．せん断速度が1/4くらいのところまで，あまり目立たないがプロットはプラトーを示している．その後は減少する．せん断速度の値が高くなると，見かけの粘度に対するの影響力は小さくなり，せん断速度がさらに大きくなると以下のプロットで示されるようにほぼ無視できるほどになっていく．

[0, 100]の区間での

のさまざまな値に対する見かけの粘度を比較するプロットを作成する：

せん断速度が大きい値になると，プロットに第二プラトーが見られる．せん断速度がゼロのときの粘度も第二プラトーに達する速さに影響を与える．

その一方でパラメータ μ_∞はせん断速度が大きいときの粘度の挙動を記述する．これは，せん断速度が無限大になる場合に見かけの粘度が収束する値である．

基準値を設定して，パラメータ値

をテストする：

[0, 3]の区間での

のさまざまな値に対する見かけの粘度を比較するプロットを作成する：

プロットから，せん断速度が増加するに従って見かけの粘度 μ_appに対する μ_∞の影響が大きくなることが分かる．せん断速度の値をもっと大きくすることでよりよい洞察が得られる．

[0, 100]の区間での

のさまざまな値に対する見かけの粘度を比較するプロットを作成する：

このプロットから，せん断速度が十分大きければ μ_∞の値が重要な役割を果たしていることが分かる．これは無限での見かけの粘度の値を制御する．

パラメータはベキ乗則の指数を表す．に近い値では，モデルはニュートン流体の動作に似ている．の値が1より小さいと，見かけの粘度はになる．の値が1より大きいと，モデルはベキ乗則モデルと類似した動作を見せ，見かけの粘度は ±∞になる．ここで符号はの符号に依存する．

基準値を設定して，パラメータ値をテストする：

[0, 3]の区間での

のさまざまな値に対する見かけの粘度を比較するプロットを作成する：

パラメータがどのくらいの速さで見かけの粘度を増加または減少させるかを制御しているのが分かる．の影響は，せん断速度が上昇するにつれて大きくなる．

[0, 100]の区間での

のさまざまな値に対する見かけの粘度を比較するプロットを作成する：

このプロットでは，見かけの粘度がに収束する速さに影響を与えているのがベキ乗則指数であることが分かる．プロット内では，に関連付けられた曲線でさえも，すべての曲線がに収束することに注意する．の値が1に近いので収束が遅いだけである．

パラメータは緩和時間である．これはモデルがに近付く速さを制御する．

基準値を設定して，パラメータ値

をテストする：

[0, 3]の区間での

のさまざまな値に対する見かけの粘度を比較するプロットを作成する：

ベキ乗則指数とは異なり、緩和時間はせん断速度が非常に小さいときでも重要な役割を果たす．

[0, 100]の区間での

のさまざまな値に対する見かけの粘度を比較するプロットを作成する：

緩和時間は見かけの粘度がに近付く速さを制御する．これはベキ乗則指数に似ているが，異なる点はが小さい値のせん断速度にも大きな影響を与えるという点である．

最後のパラメータはニュートン挙動から非ニュートン挙動への移行の速度を制御する．の値を大きくすると，モデルがニュートン挙動から非ニュートン挙動に移行するときのせん断速度の閾値が大きくなる．

基準値を設定して，パラメータ値をテストする：

[0, 3]の区間での

のさまざまな値に対する見かけの粘度を比較するプロットを作成する：

パラメータ a が第一プラトーの長さを制御しているのが分かる．しかし，一旦閾値を超えて流体が非ニュートン領域に入ると，見かけの粘度は素早く同じ値に収束する．

[0, 100]の区間での

のさまざまな値に対する見かけの粘度を比較するプロットを作成する：

せん断速度の値が大きいときの a の影響は無視できる．

クロス

クロスモデルは実験的モデルであり，ポリマー溶解や高分子溶液に便利である．クロスベキ乗則をモデル化するためにCarreauモデルを使うことができる：

ここでであり，と設定することで指定できる．

Bingham-Papanastasiou

Bingham-Papanastasiouモデルは粘塑性材料に便利である．Bingham-Papanastasiouモデルは以下を実装する：

Bingham-Papanastasiouモデルのパラメータ名：

非ニュートン流れに対する実際のBingham-Papanastasiouモデルの実装は以下に示す通りである．このモデルの構築方法は，別の非ニュートン流れモデルの構築の例として使うことができる．

Bingham-Papanastasiou非ニュートン流れモデルの実装：

ApparentViscosity["Bingham-Papanastasiou", vars_, pars_, data__]:=              
Module[                                                                         
    {model, mup, sigma, m, strainRate, shearRate, apparentViscosity},
    model = pars["FluidDynamicsMaterialModel"]["Bingham-Papanastasiou"];
    mup = model["PlasticViscosity"];
    sigma = model["YieldStress"];
    m = model["ShearRateFactor"];
    strainRate = "StrainRate" /. data;
    shearRate = PDEModels`ShearRate[strainRate + $MachineEpsilon];
    apparentViscosity = mup + sigma/shearRate * (1 - Exp[-m * shearRate]);
    apparentViscosity
]

Herschel-Bulkley-Papanastasiou

Herschel-Bulkley-Papanastasiouモデルは以下を実装する：

このモデルはベキ乗則を使ってBingham-Papanastasiouモデルの塑性粘度を計算する．ベキ乗則とBingham-Papanastasiouモデルのパラメータは，ベキ乗則で計算される"PlasticViscosity"を除いてすべて指定することができる．"MinimalShearRate"は-∞に設定される．

Herschel-Bulkley-Papanastasiouモデルのパラメータ名：

Herschel-Bulkley-Papanastasiou非ニュートン流れモデルの実装：

ApparentViscosity["Herschel-Bulkley-Papanastasiou", vars_, pars_, data__]:=
Module[
{model, powerLawPars, mup, binghamPapanastasiouPars, apparentViscosity},
    model = pars["FluidDynamicsMaterialModel"]["Herschel-Bulkley-Papanastasiou"];
    powerLawPars = pars;
    powerLawPars["FluidDynamicsMaterialModel"] = <|"PowerLaw" -> model|>;
    mup = ApparentViscosity["PowerLaw", vars, powerLawPars, data];
    binghamPapanastasiouPars = pars;
    model["PlasticViscosity"] = mup;
    binghamPapanastasiouPars["FluidDynamicsMaterialModel"] =
        <|"Bingham-Papanastasiou"-> model|>;
    apparentViscosity = ApparentViscosity["Bingham-Papanastasiou", vars, 
        binghamPapanastasiouPars, data];
    apparentViscosity
]

カスタムの見かけ粘度モデル

好きなモデルがこれまでになかった場合は，カスタムの粘度モデルを指定することで使うことができる．

カスタムの粘度モデルを指定する：

関数 myApparentViscosity は以下のシグネチャを持つ：

myApparentViscosity[_, vars_, pars_, data__]:=
Module[{},
....
]

関数 myApparentViscosity は見かけ粘度に対してスカラー値を返す．このモノグラフで紹介されている他のモデルは，見かけ粘度関数を書くためのテンプレートとして使うことができる．

カスタムの粘性応力テンソル

すべての粘性応力テンソル

は関数 myViscousStressTensor に置き換えることができる．

カスタムの粘性応力テンソルを指定する：

関数 myViscousStressTensor は以下のシグネチャを持つ：

myViscousStressTensor[vars_, pars_, data__]:=
Module[{},
....
]

関数 myViscousStressTensor は粘性応力テンソルを返す．

粘弾性流体の流れ

粘弾性流体の流れでは，非圧縮性流れに対する粘性応力テンソル（単位[Pa]）は以下で与えられる：

ここでは粘弾性応力テンソルである．粘弾性流体の流れを含む方程式はすぐに数値的に不安定になるため，安定化手法が必要になる．このバージョンのNDSolveには安定化手法が装備されていないので，粘弾性流体の流れは現在実装されない．

エネルギー輸送 - 非等温流れ

温度差がどのように流体の流れに影響を及ぼすのか，あるいは熱がどのように流体流れとともに輸送されるのかについて理解したいと思うことがあるだろう．このために，ナビエストークス方程式と熱移動方程式を結合する：

ここでは絶対温度（単位[]），は特定の圧力における非熱容量（単位[]），は熱伝導率（単位[]），は熱源密度（単位[]）である．熱方程式の詳細および導出については「熱移動」モノグラフをご覧いただきたい．熱源の中には粘性加熱または圧力変化によって記述される項がある．どちらも無視できることが多いが必要であれば加えられる．

重要な点は温度変化が密度の変化を引き起こすということである．体積が増加すると重量が減少し流体が上昇する．温度依存の浮力は上昇する．

ブシネスク近似

熱方程式とナビエストークス方程式の結合をモデル化する最も簡単なものがブシネスク近似であり，以下の仮定に基づいている[2, p.126]：

浮力項以外では密度は一定である

他の流体特性はすべて一定であると仮定する

粘性散逸は無視できる

最初の項目により，連続方程式の非圧縮形式が使える：

よって運動方程式は以下で与えられる：

が導入されており，は重力場である．積は浮力寄与である．さらに簡約するために，ブシネスク近似は密度と温度の間に線形関係を仮定する：

ここでは熱膨張係数，とは基準値である．実際，方程式はにおける一次テイラー展開である．との間の線形関係は [4]の純水のような特別な場合には二次関数に近いもので置換する必要がある．

ブシネスク近似は流れ方程式と熱移動モデルの間の結合を提供する．温度差は流体流れを駆動することができる．流体流れが等温ならばということであり，流体は静止したままである．しかし一つの壁面が熱せられると，浮力によって流体流れができる．

ブシネスク近似を使うときは注意が必要である．変化する材料パラメータは浮力項の密度だけである．つまり，その他の密度，動粘度，熱膨張係数，非熱容量，熱伝導率はすべて区間で一定なのである．ブシネスク近似は圧力は一定であると仮定する．浮力項の密度はが展開点であるテイラー展開である．

上で指摘した制約は非常に厳しいので，条件が満足していることを確認する必要がある．以下の分類を導出した[2; 5]：

ここではにおける最大温度変動であり，は代表長さである．

以下の例は水に対するこの手法を表す[2, p. 133]．近似がいつ有効であるかを推定するために，上からの不等式をコード化する．

ブシネスク近似の妥当性不等式をコード化する：

およびの水では以下のパラメータになる[5]：

規則のリストとして値を設定する：

不等式への値の挿入を簡約するためにヘルパー関数を書く：

不等式を還元する：

ブシネスク近似の水の

と

の値に対する値が有効な面積を可視化する：

圧力を考えるときも同じことができる[5]．材料定数を含む数量を積分する場合は，浮力項で使われる温度依存の密度ではなく，のような基準値を使う必要がある[4]．

レイリー・ベナール対流

ブシネスク近似の応用として，レイリー・ベナール対流を紹介する．この例では流体で満たされた矩形領域を使う．壁面は断熱されており上面と底面はそれぞれ冷却・加熱されている．まず領域といくつかのパラメータを指定する．

矩形の縦横比は4:1である．

パラメータと矩形領域を指定する：

領域を可視化する：

メッシュを生成し可視化する：

ブシネスク近似を使う熱方程式と結合した粘性ナビエストークス方程式を設定する．指定された材料パラメータを使う．

流体変数および熱移動変数を設定する：

空気の材料パラメータを指定する：

熱方程式と層流のパラメータを設定する：

ブシネスク近似を設定する：

熱方程式と組み合わせたナビエストークス方程式を設定する：

境界条件と初期条件を設定する必要がある．

境界のすべての壁面に速度に対する滑りなし境界条件を設定する：

基準圧力点を設定する：

上面と底面の間の温度差を指定する：

境界条件のパラメータを置換する：

系が静止しているように初期条件を設定する：

調整されたメッシュを使い，速度とおよび温度を二次で補間し圧力を一次で補間しながら，時間積分の進行状況とPDEを解くのにかかる合計時間を監視する．

方程式の時間積分を監視する：

最後のステップは後処理である．

従属変数に従って結果を分割する：

領域の境界の可視化を構築する：

圧力分布を可視化する：

温度分布を可視化する：

速度場を可視化する：

温度変化と速度流線をアニメーションにする：

流体流れの熱移動の例はモデルコレクションの「浮力流れ」で見られる．ここではレイノルズ数ベースのモデルが使われている．2つ目の例「熱交換器」では標準モデルが使われている．「熱分解」は熱，質量，流体の輸送の組合せを示している．

流れ境界条件

流体流れモデルの幾何学的領域は有限である．領域の境界では流れが埋め込まれている環境をモデル化する境界条件が必要である．境界条件は境界の壁面や領域への流入，領域からの流出を記述する．境界条件を簡単に説明するためには，いくつかの用語を導入するのが便利である．が境界に垂直な速度成分を表し，が境界に接線方向である速度成分を表すとしよう．とは法線方向における両方の成分の導関数を表す．

先に紹介したリッドドリブンキャビティの例のように，軸と軸に平行な2D形状を考える．ここでは垂直境界と水平境界の2つがある．垂直境界では：

であり，下が成り立つ：

ピンクで示された垂直境界では，垂直流体流れ成分と水平流体流れ成分がある．

水平境界では：

であり，以下が成り立つ：

ピンクで示された水平境界では，垂直流体流れ成分と水平流体流れ成分がある．

座標軸に平行ではない境界の場合は，成分とをとから計算する必要がある．

流入条件

流入とは入口とも言われる．両方の速度成分が規定される：

ここでとは速度ベクトルに垂直，水平な成分である．流入条件を指定するためにとが与えられる．通常の流入プロファイルは以下である：

ここでは垂直方向の流体速度を指定する．

ここで覚えておかなければならないのは，流入条件は隣接する境界条件と矛盾していてはならないということである．例えば，入口の境界における流入プロファイルが0ならば，滑りなし壁面境界条件を使うことができる．入口の境界における流入プロファイルが非零ならば，隣接する滑りなし壁面境界条件を使うと，ソルバがモデルを収束するのが不可能ではないとしても難しくなる．場合によっては境界の流入条件が非零のときは，滑り壁面境界条件を使った方がよいことがある．

Fitを使うと放物型流入プロファイルが簡単に作成できる．の流入と最大流入速度1[]が与えられたとする．

端点での速度0，中央の速度1[

]のデータ点を作成する：

放物型流入プロファイルを作成する：

放物型流入プロファイルを可視化する：

流出条件

流出は出口とも呼ばれる．両方の速度成分は法線に対して一定のままである：

壁面条件

壁面は流体を含む任意の不透水境界である．壁面境界での条件の一つは，流体は壁面を通過できないということである．したがって，最小必要条件は以下になる：

この最小必要条件が満たされたときだけ滑り条件と言う．2つ目の条件は，追加的に壁面条件に接線方向の流体流れはないというものである．この2つ目の場合は滑りなし条件と言う．

壁面境界条件が流入境界条件と一致していることに注意を払う必要がある．詳細は「流入条件」のセクションを参照のこと．

滑りなし

通常の壁面境界条件は滑りなし条件である．滑りなしとはまず壁面を通過する流体がなく，次に流体が壁面に付着しているというものである．2つ目の条件は境界の摩擦をモデル化する．つまり以下を意味する：

よって，流体流れ速度ベクトルはこの境界では0に設定される．

これはDirichletConditionで設定することができる．

滑りなし壁面境界を指定する：

滑りなしの例

この例は滑りなし壁面境界条件の応用を示す．2つの半円板が切り取られた矩形の流れ領域がある．左側は流入プロファイルがあり右側は流出条件がある．領域の残りの部分は滑りなし境界条件としてモデル化された壁面である．

領域を設定し，メッシュ化して可視化する：

層流方程式を設定する：

流入および流出の条件を指定する：

滑りなし壁面境界条件を指定する：

境界条件を集める：

流体流れPDEを解く：

流れ場を可視化する：

における

速度成分をプロットする：

およびにおいて速度は0である．

および

における

速度成分をプロットする：

およびにおいて速度は0である．

滑り

滑りにおいても壁面を通過する流体はないが，滑りなし条件と異なり摩擦による損失はない．つまり以下が言える：

座標軸に平行な場合，垂直壁では以下となる：

座標軸に平行な垂直壁面について滑り壁面境界を指定する：

軸に垂直な場合，水平壁では以下となる：

座標軸に平行な平行壁面について滑り壁面境界を指定する：

牽引

ナビエストークスの運動方程式を考えてみよう：

粘性応力テンソルは以下で与えられる：

ここでひずみ速度は以下である：

すべての項を展開すると次が得られる：

ここでは全応力である．全応力とは粘性応力に圧力成分を足したものである．

牽引を規定するために次を指定する：

ここでは牽引値ベクトルである．すべての項を置き換えると以下が得られる：

成分形式の全応力テンソルは次である：

境界条件のタイプはNeumannValueで指定することができる．

牽引の例

次の例はアニュラス領域における規定された牽引境界条件を持つ流れを考える．辺の外側には滑りなし境界条件がある．これにより方向と方向の流体速度が0，に設定される．

内側の境界では，速度を定義する代りに牽引を規定する．応力ベクトルの牽引は次で与えられる：

ここではとして与えられる単位半径ベクトル，はで与えられる単位接線ベクトルである．とはそれぞれ方向方向の直交単位ベクトルである．この例ではの法線方向との接線方向に牽引を設定する．つまり，内側境界に接線方向の流れだけで，横切る流れはない．

パラメータを設定する：

形状と調整されたメッシュを指定する：

層流演算子を指定数する：

外側境界には滑りなし条件がある．

外側境界に滑りなし条件を設定する：

内面ではクロス積を持つ外向きの単位法線から接線単位法線を計算することができる．

クロス積を使って単位法線ベクトル

から単位接線を計算する：

方向，方向のクロス積の第1成分と第2成分をそれぞれ抽出するためにIndexedを使う．

牽引境界条件を設定する：

同様に，内面の単位法線ベクトルはであり単位接線はなので，接線単位法線はとして与えることもできる．

牽引境界条件と同等の集合：

圧力についての境界条件を指定しないと，圧力値は浮き，NDSolveは十分な境界条件が指定されていないという警告を発する．

方程式を解く：

圧力の解と速度場を可視化する：

牽引境界条件が動作することを検証するためには，速度の動径成分と接線成分を計算し，それをすべてのについてでプロットすることができる．動径速度についてはに比例する解が想定され，接線速度についてはに比例する解が想定される．

速度の動径成分と接線成分を計算する関数を定義する：

に比例する速度の動径成分と関数をプロットする：

に比例する速度の接線成分と関数をプロットする：

次に内側境界における法線応力とせん断応力を計算し，それらが規定の牽引境界条件と一致するかをチェックする．全応力テンソルは以下で与えられる：

牽引ベクトルの法線および接線成分は次で計算できる：

および

応力テンソル

を計算する関数を作成する：

単位法線ベクトルを計算する：

接線ベクトルのための関数を定義する：

法線応力を計算する関数を定義する：

せん断応力を計算する関数を定義する：

計算された法線応力とせん断応力を

で可視化する：

における法線応力・せん断応力と境界条件との差分を可視化する：

初期条件

時間依存問題では初期条件を指定する必要がある．

初期条件は連続方程式を満足しなければならない．

流体流れモデルの収束

流体流れPDEは収束が難しいことがある．このセクションはNDSolveが流体流れモデルの解を見付けられないときに試すヒントを提供する．ここで示すメソッドを組み合せることもできる．

初期シード（種）

通常非線形のソルバには最初の推測，つまり開始するためのシードが必要である．このシードが解から遠すぎるとソルバは収束しないことがある．ソルバが収束する半径のことを収束半径と言う．初期シードが近いほど最終的な解はよくなる．

反復ステップ

ここでは小さいレイノルズ数の流れ場の解を使って，大きいレイノルズ数の計算のシードにする．このためにパラメトリック関数を使う．

リッドドリブンキャビティの例を設定する：

メッシュを設定する：

従属変数を指定することでソルバが最適化された方程式系を作成することができる．

DependentVariablesを持つ構築されたメッシュを使って，新しいパラメトリック関数を作成する：

初期の解0で20ステップでレイノルズ数3万までの解を計算する：

大きいレイノルズ数の流れ場について反復的に解く：

レイノルズ数3万における解のベクトルプロット：

方程式の変更

最初はばかばかしく聞こえるかもしれないが，なぜ方程式を変更するのだろうか？より安定したより還元された方程式集合を使うことが可能な場合がある．それからその解を完全な方程式のシードとして使うことができる．

ストークス方程式

初期シードが指定されていないと，ソルバはすべての従属変数について0を仮定する．流体力学の場合，速度場に対して初期速度0が使われ圧力が使われると，実際ストークスタイプの方程式を解いていることになる．

参考文献

1. Simple Fields of Physics; Backstrom, G.; 2006; GB Publishing; ISBN 91-975553-0-4

2. Numerische Simulation in der Strömungsmechanik; Griebel M. et al.; 1995; Vieweg; ISBN 3-528-06761-6

3. The Efficient Engineer; https://www.youtube.com/watch?v=VvDJyhYSJv8; retrieved 2022/12/05

4. On the use and misuse of the Oberbeck–Boussinesq approximation; Barletta, A.; Celli, M; Rees, D.A.S; https://arxiv.org/abs/2202.10981v2; retrieved 2023/03/21

5. The validity of the Boussinesq approximation for liquids and gases; Gray, D.; Giorini, A.; International Journal of Heat and Mass Transfer Volume 19, Issue 5, May 1976, Pages 545-551; https://doi.org/10.1016/0017-9310(76)90168-X

6. “High-Re Solutions for Incompressible Flow Using the Navier-Stokes Equations and a Multigrid Method; Ghia, U., Ghia, K.N., Shin, C.T.; Journal of Computational Physics vol. 48, pp. 387–411, 1982.; https://doi.org/10.1016/0021-9991(82)90058-4

トップへ

層流

目次

はじめに

概要的な例と解析タイプ

定常解析

後処理

一般的な可視化テクニック

パラメトリック解析

時間依存解析

方程式

ニュートン流れと非ニュートン流れ

ベキ乗則

ベキ乗則の例

Carreau

クロス

Bingham-Papanastasiou

Herschel-Bulkley-Papanastasiou

カスタムの見かけ粘度モデル

カスタムの粘性応力テンソル

粘弾性流体の流れ

エネルギー輸送 - 非等温流れ

ブシネスク近似

レイリー・ベナール対流

流れ境界条件

流入条件

流出条件

壁面条件

滑りなし

滑りなしの例

滑り

牽引

牽引の例

初期条件

流体流れモデルの収束

初期シード（種）

反復ステップ

方程式の変更

ストークス方程式

参考文献

関連するテクニカルノート