分子動力学計算による複雑流体中のカルマン渦の解析

[4]

 最近の研究から「分子動力学計算による複雑流体中のカルマン渦の解析」

東京大学 物性研究所 浅野優太 [email protected] 慶應義塾大学 理工学部 渡辺宙志 [email protected] 東京大学 物性研究所 野口博司 [email protected] 概要  流体に極微量の高分子を加えることで流動抵抗が劇的に低減するトムズ効果や，水中で回転するタービンの周 りでキャビテーションにより発生する気泡の運動など，工学上重要な問題はミクロスケールからマクロスケール までシームレスに影響を与えるものが多い．これらの問題は典型的なマルチスケール・マルチフィジクス問題で あり，解析が極めて難しく，メカニズムはよく分かっていない．計算機の発達により，流れそのものを分子動力 学 (MD) 法を用いた全粒子計算で解析するという力任せな方法が現実味を帯びてきた．本稿では，複雑流体の流 動，キャビテーションを伴う流動に関して，著者らが行ってきた MD 計算の紹介を行う． キーワード：分子動力学法，流体，カルマン渦，高分子，キャビテーション 1 はじめに 流体の流れは，エアコンから吹き出る空気や水道管 内の流れのような身近なものから，エネルギー輸送を 担うパイプラインの流れやタービンなどのエネルギー 供給システムのような工業プロセスに至るまで，幅広 い分野で非常に重要な役割を果たしている．したがっ て，流れの特性を良く理解し，流れを適切に制御する ことは，流体力学・流体工学において重要な研究課題 である．特に，工業的な流れの制御については，高分 子や界面活性剤などを含む複雑流体の性質を利用した 流動抵抗の低減メカニズムや，キャビテーションなど の相転移を伴う気液二相流による熱・物質移動といっ た，複雑な流動現象の解明が望まれている．そのため， これまでに多数の実験や Navier-Stokes 方程式に基づ く計算機シミュレーションなど，様々なアプローチに よる解析が行われてきた．しかしながら，複雑流体や 気液二相流などの工学上重要な問題は，ミクロスケー ルのダイナミクスとマクロスケールの流動が共存する マルチスケール・マルチフィジクス問題であり，その 解析は極めて難しい． このような困難を克服し，ミクロスケールのダイナミ クスを捉える新しい手法として，格子ボルツマン法 [1] や，multi-particle collision dynamics [2] といったメソ スケールの解析や，マクロスケールの流動とミクロス ケールのダイナミクスの連成解析 [3] などが開発されて きた．しかしながら，今日の計算機能力を考慮すれば， 分子動力学 (MD) 法による全粒子計算によって，流れ 場そのものを取り扱うことが可能であると考えられる． MD 計算による流体解析では，基本的には，与えるも のが粒子間相互作用だけである．そのため，相転移や 高分子のコンフォメーションの変化など，他の手法で は状態方程式や構成方程式などの現象論的なモデルが 必要であった現象が自然な形で発生する．その代わり， 流動の特徴を解像できる規模の計算セルを用意しなけ ればならないため，必然的に規模の大きな計算になる． 著者らはこれまで，MD 計算によって，高分子の添加 や相転移による流れ場への影響を分子スケールから調 べてきた．本稿では，これまでの取り組みについて簡 潔に紹介する． 以下の 2 節では，粒子間相互作用に共通の関数形を 用いているため，予め示しておく． u(r) = { u0(r)− u0(rc) + (r− rc)u′0(rc) (r≤ rc) 0 (r > rc) (1) u0(r) = 4ϵ [(_σ r )12 −(σ r )6] (2) ここで，r は粒子間距離であり，ϵ と σ は，それぞれ， エネルギーと長さのスケールを表す．プライムは r に ついての微分を表す．rcはポテンシャル関数の切断距 離である．以下では，全ての物理量を，エネルギーを ϵ，長さを σ，時間を τ = σ√m/ϵ の単位でそれぞれ測 る．m は粒子の質量である．

2 高分子添加によるカルマン渦の抑制 2.1 背景 ニュートン流体に極微量の高分子を添加すると，流 動抵抗が劇的に低減するトムズ効果 [4] や，カルマン 渦の放出周期が長くなること [5] など，流体の性質が 大きく変わることが実験的に知られている．多くの工 業的な流れでは，渦の放出や，乱流状態になっている ことが多く，高分子の添加による流れ場の制御が注目 されている．しかしながら，流れ場の中での高分子の レオロジーが十分に理解されていないため，高分子の 添加によって生じる現象のメカニズムはよくわかって いない． 高分子の添加によって生じる流動特性の変化は，流 れ場で発生している渦と高分子の相互作用が重要な役 割を果たしていると考えられている [6]．そこで，渦を 伴う流れとして代表的であるとともに，周期性から解 析に適しているカルマン渦を対象に選んだ．渦と高分 子の関係を分子スケールから議論するために，高分子 の添加によるカルマン渦への影響を MD 計算によって 解析した [7]． 2.2 計算方法 実験 [8] を参考にして，2 次元の MD 計算を行い， ニュートン流体と希薄高分子溶液の円柱周りの流れの解 析を行った．溶媒粒子間の相互作用には，(1) 式において rc = 21/6とした，Weeks-Chandler-Andersen(WCA) 相互作用 [9] を用いた．ここでは，WCA 粒子のみで構成 される溶液をニュートン流体と呼ぶ．高分子は Kremer-Grest モデル [10] を用いて表現した．図 1 に示すよう に，計算セルは大きさが Lx× Ly = 1000× 500 の長方 形とし，全方向に周期境界条件を課した．図中の黒丸 は円柱を表す．円柱は WCA 粒子をその表面上に固定 することで表現した．円柱の直径は D = 100 であり， 中心の位置は (x, y) = (250, 250) とした．x 軸方向に一 様な流入速度 V を与えるために，系の一部に Langevin 熱浴によって流速を制御する領域を設けた (図 1 の斜 線部)．この熱浴は円柱の後ろに発生したカルマン渦な どを消す役割も担っている．温度は T = 1 とし，流入 速度 V は音速 vs ≃ 6 を越えない範囲 (1.5 ≤ V ≤ 3) で変化させた．溶液の密度は ρ = 0.83 とし，高分子は 短鎖長 (セグメント数 Ns = 10) と長鎖長 (Ns = 100) の 2 種類を用いた．高分子の濃度は絡み合い濃度以下 とし，最大濃度は体積分率で ϕ = 0.107 とした．全粒 子数は N = 408 482 である．流れ場を特徴づける無次 元量のレイノルズ数 Re は，ニュートン流体の粘性係 数 η を用いて，Re = ρDV /η と評価した． 0 250 500 0 250 500 750 1000 y x Cylinder 図 1: 計算セルの概略図．図中の黒丸は直径 D = 100 の円柱であり，中心位置は (x, y) = (250, 250) である． 斜線部 (800≤ x ≤ 1000) ではランジュバン熱浴により 温度と流速を制御している． 2.3 結果と考察 図 2 に渦度場 ω = ∂vy ∂x − ∂vx ∂y のスナップショットを 示す．vxと vyは，それぞれ，流速の x 成分と y 成分 である．ここでは計算セルを 10× 10 の格子に分割し， 中心差分によって各格子の ω を求めた．ニュートン流 体では (図 2(a))，円柱の背後にカルマン渦が発生して いることがわかる．一方，長鎖長の高分子を添加する と，渦がぼやけることがわかる (図 2(b))．短鎖長の高 分子の場合は，渦度場の様子はニュートン流体と変わ らず，粘性を変える効果しか現れなかった．したがっ て，高分子の鎖長が重要な役割を果たしていることが わかる． カルマン渦が発生すると，円柱には渦の放出周期の 揚力 FLが働くため，そのフーリエ変換から渦の放出周 期が得られる．図 3 に揚力係数 CL= 2FL/(ρDLzV2) のフーリエ変換を示す．Lzは z 軸方向の厚みであり， 2 次元系の場合は便宜上 Lz= 1 とする．図 3(a) に示す ように，ニュートン流体の場合，渦の放出周波数の位 置に鋭いピークがあり，周期的に渦の放出が起こって いることがわかる．短鎖長の高分子を添加した場合も， ニュートン流体と同様に鋭いピークを示していた (図 3(b))．ところが，長鎖長の高分子を添加した場合は， 高分子濃度の増加に伴って，ピーク位置が低周波数側 にシフトし，スペクトルの幅が広がることがわかった (図 3(c) と図 3(d))．MD 計算で得られた，(1) 渦のぼ やけ，(2) 周波数の低下，(3) 鎖長の影響は，実験で知 られている結果 [8] とよく一致する． 流れ場の中での高分子の振る舞いを調べるために， 各位置での慣性半径 Rgを評価した．慣性半径は次式 で与えられる： R2_g = 1 Ns Ns ∑ i=1 ⟨( ⃗ Ri− ⃗RG )2⟩ . (3)

100 250 400 550 700 x 0 250 500 y -0.1 -0.05 0 0.05 0.1 Vorticity

(a) Newtonian fluid

100 250 400 550 700 x 0 250 500 y -0.1 -0.05 0 0.05 0.1 Vorticity

(b) Dilute polymer solution

図 2: (a) ニュートン流体と (b) セグメント数 Ns= 100 の高分子を添加した溶液の渦度場のスナップショット． レイノルズ数は Re = 64，高分子の体積分率は ϕ = 0.107 である． ここで， ⃗Riと ⃗RGは，それぞれ，高分子の i 番目のセ グメントと重心の位置ベクトルである．_{⟨· · · ⟩ は統計平} 均を表す．静止流体中の短鎖長と長鎖長の高分子の慣 性半径は，それぞれ，Rg= 1.8 と 10.2 である． 図 4(a) に示すように，短鎖長の高分子の慣性半径は ほぼ静止流体中の値となっている．したがって，高分 子は円柱通過時に受ける強いせん断や，円柱後方の渦 の中で変形せず，粒子のように振る舞っていることが わかる．一方，長鎖長の高分子は，円柱近傍，及び円 柱後方の渦の中で大きく引き伸ばされていることがわ かる．また，レイノルズ数 Re = 50 では円柱の後方で 僅かに減少傾向を示しているが，Re = 64 ではほぼ一 定値である．Re が大きい場合，発生する渦の強さも大 きくなるため，渦列の付近の高分子を巻き込んでいる ためであると考えられる．高分子の運動と渦の強さの 関係を調べるために，各位置での高分子の配向秩序 Q を求めた．配向秩序は次式によって計算した： Q = 2⟨cos2θ⟩− 1. (4) ここで，θ は末端間ベクトルと x 軸の成す角である．高 分子が流れ方向に配向している場合は Q = 1，流れに 垂直方向に配向している場合は Q =−1，ランダムな 0 0.1 0.2 0.3 AL

(a) Newtonian fluid

0 0.1 0.2 AL (b) N_s=10 φ=0.107 0 0.1 0.2 AL (c) N_s=100 φ=0.024 0 0.1 0.2 0.15 0.2 0.25 0.3 AL f ~ (d) Ns=100 φ=0.107 図 3: レイノルズ数 Re = 64 における，円柱に働く揚 力係数の振幅 AL：(a) ニュートン流体，(b) 高分子溶液 (高分子のセグメント数 Ns= 10，体積分率 ϕ = 0.107)， (c) 高分子溶液 (Ns= 100，ϕ = 0.024)，(d) 高分子溶液 (Ns= 100，ϕ = 0.107)．実線はガウス関数によりフィッ トしたものである．破線はニュートン流体のピーク位 置を表す．横軸は無次元化した周波数 ˜f = f D/V (f ： 揚力係数の周波数，D：円柱直径，V：流入流速) である． 場合は Q = 0 となる．図 4(b) に配向秩序の各位置で の値を示す．短鎖長の高分子の場合，ほぼ一様な値で あり，配向もほぼランダムである．しかしながら，長 鎖長の高分子の場合，不均一な値となっており，円柱 近傍，及び円柱後方の渦列中で特徴が現れている．円 柱通過時に引き伸ばされ，流れ方向に配向した高分子 が，円柱の後方の渦の中でも流れ方向に配向する傾向 にあることがわかる．円柱から離れていくにしたがっ て，配向秩序は失われていくが，レイノルズ数が高い 方がその程度が大きい．これは渦によって，周囲の高 分子が巻き込まれているためである． せん断による高分子の伸長と渦による高分子の巻き 込みが，渦の抑制に重要な役割を果たしていることを 分子スケールから明かにした．また，本稿では詳細に は立ち入らないが，カルマン渦における有限サイズ効 果を調べ，計算領域のサイズと粘性の密度依存性によ る高密度流体での粘性の非一様性の影響を明らかにし， ニュートン流体の振る舞いが再現できる条件を明確に

した．これについては文献 [11] を参照されたい． 0.8 1.2 1.6 Rg / R eq g N_s = 100, Re = 64 N_s = 100, Re = 50 (a) -0.3 0 0.3 150 250 350 450 550 650 750 Q x N_s = 10 , Re = 64 N_s = 10 , Re = 50 (b) 図 4: レイノルズ数 Re = 50 と 64 における，短鎖長 (高分子セグメント数 Ns= 10) と長鎖長 (Ns= 100) 高 分子の (a) 慣性半径 Rgと (b) 配向秩序 Q の各位置 x での平均値．Req g は静止流体中での高分子の慣性半径 である．各量は時間と y 軸方向で平均化してある．高 分子の体積分率は ϕ = 0.107 である． 3 カルマン渦へのキャビテーションの影響 3.1 背景 キャビテーションは，高流速の液体中の局所的な圧 力変化によって生じる，気泡の生成・成長・消滅を伴 う流動現象である [12]．船舶の推進機の性能低下，騒 音や振動，エロージョンなど，キャビテーションは流 体機械に大きな影響を与えるため，その発生メカニズ ムを理解することは工学的に非常に重要である．しか しながら，マクロな流動場中でのミクロな気泡核のダ イナミクスを議論することが困難なため，キャビテー ションの発生メカニズムは未だ十分な理解がなされて いない．本研究では，円柱周りの流れで発生するキャ ビテーションの気泡の生成成長過程，及び周囲の流れ に与える影響を MD 計算によって調べた [13]． 3.2 計算方法 3 次元の MD 計算により，相転移を含むニュートン流 体の円柱周りの流れの解析を行った．粒子間相互作用は， (1) 式において，rc= 2.5 とした Lennard-Jones(LJ) 相 互作用を用いた．計算セルは大きさが Lx× Ly× Lz= 3600×1000×100 の直方体とした (図 5)．全方向に周期 境界条件を課した．2 つの円柱を並べて周期境界条件を 課すことで，上下の渦列が非対称な構造も取り得るよう になる．前節と同様に，系の一部には流速制御，平衡化 のための熱浴領域を設けた (図 5 の斜線部)．流入流速は V = 1 とした．密度は ρ = 0.4 とし，温度を 1≤ T ≤ 2 の範囲で変化させた．全粒子数は N = 143 421 343 で ある．臨界点の温度と密度は，それぞれ，Tc = 0.94， ρc= 0.32 であり，降温により臨界点に近づけていくと 円柱近傍での状態変化によって気泡が生成される． 0 250 500 750 1000 0 1000 2000 3000 y x Cylinder 1 Cylinder 2 3600 図 5: 計算セルの概略図．直径 D = 100 の円柱 1 と 2 は，それぞれ，(x, y) = (1000, 250) と (1000, 750) に配 置した．斜線部 (3200≤ x ≤ 3600) ではランジュバン 熱浴により温度と流速を制御している． 相転移の有無による違いを議論するために，WCA 流体による単相流の解析も行った．また，円柱間の渦 列の位相による影響を調べるために，1 つの円柱のみ を置いた場合の解析も行った．その場合は図 5 の上下 のどちらか半分の領域を消去した計算セルになる． 3.3 結果と考察 図 6 に密度場のスナップショットを示す．高温 (T = 2) では，気泡は発生せず，非キャビテーション流れになっ ている．円柱の後方には，逆位相に同期したカルマン 渦が発生した．円柱近傍，及び渦の中心部は他の領域 に比べて密度が低くなっていることが分かる．これら の領域は他に比べて圧力が低くなっているため，温度 を下げていくとキャビテーションが発生することが期 待される． 温度を下げていき，密度場の変化を調べた．T = 1.3 では，円柱近傍からカルマン渦の放出に連動して，気 泡が発生していることがわかる．局所平衡を仮定した 場合，円柱近傍の低密度領域は，気相であることがわ かった．一方，円柱後方の渦の中心部は，低密度流体 であることもわかった．しかしながら，現段階では，非 平衡状態において気液の判断をする術が無いため，密 度が低い部分を気体であると考えるのがもっともらし いと思われる．気泡発生直後のこの温度では，T = 2 と同様に逆位相同期したカルマン渦が見られ，渦の構 造に対する影響は見られない． T = 1.25 では，円柱の背後に気相領域が形成され 円柱に付着していることがわかる．さらに，渦が形成

図 6: LJ 流体の密度場のスナップショット：(a) 非キャ ビテーション流れ，(b)-(d) キャビテーション流れ． されるまでの距離が増加していることがわかる．また， 僅かではあるが，上下のカルマン渦の位相に乱れが生 じている．温度をさらに下げ，T = 1.2 になると，円柱 背後の気相領域がさらに拡大する．また，上下の渦構 造が非対称になる．この非対称構造は，長周期で入れ 替わる．WCA 流体のように，相転移が起こらない流 体の場合は，温度変化によってこのような流動場の変 化は起こらない．気泡の発生の伴う密度揺らぎによっ て，渦構造に大きく影響を及ぼすことがわかった． 渦構造が変化すると，渦の放出によって励起される 振動にも影響が現れるはずである．そこで，円柱に働 く揚力係数の振動振幅 ALと，渦形成長さ Lfのレイノ ルズ数 Re 依存性を調べた．図 7 に ALと Lfの Re 依 存性を示す．Lfは，円柱の中心軸と円柱後方において 流速の揺らぎが最大となる位置の間の距離として求め た．また，LJ 流体と WCA 流体の粘性係数は降温に 伴って減少していくため，Re は増加する．WCA 流体 の場合は Re の増加に伴って，ALは増加傾向，Lf は 減少傾向であり，ともに単調な変化を示す．ところが， LJ 流体の場合，Re に対して非単調な変化を示してい 0 0.2 0.4 0.6 0.8 AL WCA Single d=5D WCA Double d=5D Single d=10D Single d=5D Double d=5D (a) 0 500 1000 1500 90 100 110 120 Lf Re (b) 図 7: 揚力係数の振幅 ALと，渦形成長さ Lfのレイノ ルズ数 Re 依存性．実線は LJ 流体，破線は WCA 流体 の結果を示す． ることがわかる．これは，Re の増加と気泡発生の競合 によるものと考えられる．円柱間の距離 d = 5D では， 気泡は Re≃ 105 で発生する．気泡が発生すると，ま ず ALの減少が始まり，Re≃ 108 で消失する．ALの 消失と同時に，Lfが急激に増加する．d = 10D とした 場合，流路幅の増加によって円柱通過時に生じる圧力 低下が抑えられるため，気泡の発生温度が低温側 (高 Re 側) にシフトする．しかしながら，気泡発生後の振 る舞いは，d = 5D の場合と同様である． WCA 流体の場合，1 円柱と 2 円柱の場合で，ALと Lfはほぼ同じ傾向を示しているが，LJ 流体の場合は 気泡が発生する前 (95 < Re < 105) の振る舞いが異 なっている．ALには上向きのピーク，Lfには下向き のピークがあり，その位置は 2 円柱では Re≃ 101 で あり，1 円柱では Re≃ 105 である．このピークは，気 泡の発生直前に生じる非常に小さな気泡生成，或いは， 臨界点近傍であることに起因する密度揺らぎによって 生じたものと考えられる．したがって，位相の引き込 みによる振動の増幅によって，2 円柱のピーク位置は 1 円柱よりも高温側 (低 Re 側) になる．また，1 円柱 のピークが小さいのは，気泡の発生による影響によっ て打ち消されたためと考えられる．気泡発生後は，気 泡の影響により渦列間の干渉による影響が弱まるため， 位相の違いによらず，ほぼ同じ振る舞いとなる．円柱 近傍で発生した気泡により，揚力や流動場の特性が大 きく変化することを分子スケールから明らかにした．

4 おわりに MD 計算によって，高分子添加によるカルマン渦の 抑制，及びカルマン渦に対するキャビテーションの影 響を調べた．高分子添加については，長鎖長の高分子 を添加すると，せん断による高分子の伸長，及び渦に よる高分子の巻き込みが渦の抑制に重要な役割を果た していることを明かにした．キャビテーションについ ては，気泡の発生によって，渦構造が変化することや， 揚力の振動振幅が減少し，最終的に消失することを明 かにした． 従来の手法と比較すると，MD 計算による流体解析 では，構成方程式や状態方程式などの現象論的なモデ ルを使わずに現象を議論できるという点が強みである． したがって，計算コストの問題をクリアできれば，プ ロペラや翼周りの流れや，複雑流体中の気液二相流と いった様々な流れの解析に有望なツールであると考え られる．今後は，高分子と気泡発生の両者が共存する 流動場の解析に取り組み，より複雑な流動特性を分子 スケールから解明していきたい． 謝辞 本研究の遂行にあたって，東北大学の川勝年洋教授， 森井洋平氏，九州大学の津田伸一准教授，國嶋雄一氏， 東京大学物性研究所の樋口祐次氏に多数の助言を頂き ました．ここに記して，感謝の意を申し上げます．本 研究は，文部科学省ポスト「京」萌芽的課題「基礎科 学の挑戦−複合・マルチスケール問題を通した極限の 探求」，JSPJ 科研費 JP15K05201 の助成を受けたも のです．本研究の計算には，東京大学物性研究所，自 然科学研究機構計算科学研究センターのスーパーコン ピュータを利用させて頂きました． 参考文献

[1] S. Succi, The Lattice Boltzmann Equation: For Fluid Dynamics and Beyond, Oxford University Press Inc., New York (2001).

[2] A. Malevanets and R. Kapral, J. Chem. Phys.,

110, 8605 (1999).

[3] X. Xu and P. Yu, J. Non-Newtonian Fluid Mech.,

229, 27 (2016).

[4] G. E. Gadd, Nature, 206, 463 (1965). [5] G. E. Gadd, Nature, 212, 874 (1966).

[6] T. Min, J. Y. Yoo, H. Choi, and D. D. Joseph, J. Fluid Mech., 486, 213 (2003).

[7] Y. Asano, H. Watanabe, and H. Noguchi, J. Chem. Phys., 148, 144901 (2018).

[8] R. Cressman, Q. Bailey, and W. I. Goldburg, Phys. Fluids, 13, 867 (2001).

[9] J. D. Weeks, D. Chandler, and H. C. Andersen, J. Chem. Phys., 54, 5237 (1971).

[10] K. Kremer and G. S. Grest, J. Chem. Phys. 92, 5057 (1990).

[11] Y. Asano, H. Watanabe, and H. Noguchi, J. Phys. Soc. Jpn., 88, 075003 (2019).

[12] C. E. Brennen, Cavitation and Bubble Dynamics, Oxford University Press Inc., New York (1995). [13] Y. Asano, H. Watanabe, and H. Noguchi, J.

Chem. Phys., 152, 034501 (2020). 著者紹介 浅野優太 (博士 (理学))：〔経歴〕 2014 年愛媛大学大学院理工学研 究科博士課程修了，理化学研究所 計算科学研究機構，名古屋大学工 学研究科を経て,2016 年から現所 属．〔専門〕分子シミュレーション， 統計物理学．〔趣味〕音楽． 渡辺宙志 (博士 (工学))：〔経歴〕 2004 年東京大学工学研究科博士 課程修了，名古屋大学大学院情報 科学研究科，東京大学情報基盤セ ンター，東京大学物性研究所を経 て，2019 年より現職．〔専門〕統 計力学，分子動力学法．〔趣味〕プ ログラム． 野口博司 (博士 (学術))：〔経歴〕 2000 年名古屋大学人間情報学研究 科博士課程修了, 分子科学研究所, ユーリッヒ研究所を経て, 2008 年 から現所属.〔専門〕ソフトマター 物理, 分子シミュレーション．