技 術 論 文

技 術 論 文MPS 法による粘弾性流体の数値シミュレーション

MPS 法による粘弾性流体の数値シミュレーション
Numerical Simulation of Viscoelastic Fluid by the MPS Method

博士（工学）

福澤 洋平 *

Dr. Yohei Fukuzawa

博士（工学）

松永 拓也 **

Dr. Takuya Matsunaga

要

工学博士

越塚 誠一 ***

Dr. Seiichi Koshizuka

旨

粘弾性流体の挙動と自由表面形状を正確に予測することを目的とし、移動サーフェスメッシュ LSMPS 法に粘弾性流体の Oldroyd-B モデルを適用し、粘弾性流体の挙動を高精度に予測する演算手法を新たに開発した。本解析手法を用いることで、粘弾性流体の特異挙動であるワイゼンベルグ効果を高精度に再現・予測でき、開発した解析手法の有用性が示唆された。

Synopsis
We have developed a method to predict phenomena of viscoelastic fluids by particle method. With the aim of accurately
simulating the behaviors of viscoelastic fluids with free surfaces, the moving surface mesh-incorporated LSMPS（LeastSquares Moving Particle Semi-implicit）method has been extended to viscoelastic fluids based on the Oldroyd-B model. The
developed method is verified using Weissenberg effect, which is a phenomenon unique to viscoelastic fluids. As a result, wesucceeded in reproducing the Weissenberg effect.

１.

はじめに値振動が発生しやすいことから、これらを抑制して高精度に予測するための課題を有する。また、自由表面をスクリュ回転による押出機内部の流れ場はせん断流れ扱う主な手法として界面捕獲法である VOF（Volume ofが支配的となるため、流動解析を実施する際には樹脂のFluid）法やレベルセット法が用いられるが、これらの手流体モデルを一般的に純粘性非ニュートン流体として扱法は自由表面の応力境界条件を厳密に取り扱うことがでう。これらの流動解析は実績も豊富でありほぼ確立されきないため、粘弾性の効果によって大きく変形する自由た技術となっている。一方で、押出後の自由表面を有す表面形状の予測が難しい。一方で、ラグランジュ法である流動解析では、ダイスウェル現象やワイゼンベルグ効る粒子法は連続体を粒子の集合体で表し、各粒子の運動果に代表されるように、粘弾性流体の特性が顕著に現れを計算することで連続体の挙動を予測する手法であり、る。これらの挙動を正確に再現・予測するためには、流自由表面の追跡や大変形問題を容易に予測できることが体モデルを粘弾性流体として扱い、自由表面形状を妥当最大の特徴である。粒子法の解法として MPS（Movingに予測する解析手法が必要となる。粘弾性流体の流動解（Smoothed Particle
Particle Simulation）法（1）と SPH

析では格子法による予測技術が古くから開発されていHydrodynamics）法（2）が広く普及しており、粒子法の特る。格子法はオイラー法であるため支配方程式に移流項徴を活かして幅広い産業分野で実用されている。しかが存在し、この移流項を離散化する過程で数値拡散や数しながら、計算過程で粒子配置が乱雑になった際には、*：イノベーションマネジメント本部 先端技術研究所 Advanced Technology Laboratory, Innovation Management Headquarters**：東京大学大学院 新領域創成科学研究科 人間環境学専攻 Human and Engineered Environmental Studies, The University of Tokyo***：東京大学大学院 工学系研究科 システム創生学専攻 Department of Systems Innovation School of Engineering, The University of Tokyo（1）

技 術 論 文

MPS 法による粘弾性流体の数値シミュレーション

圧力振動が生じ予測精度が著しく低下する課題がある。

ここで、気体密度は無視できるほど微小と仮定すると、この課題に対して、圧力振動を低減し予測精度を高めるための研究が活発に行われるようになり、中でも最小二（一定）と（一定）となり、以下の式となる。

乗法を用いる高精度スキームの LSMPS（Least Squares（5）

Moving Particle Semi-implicit）法（3）によって予測精度が劇的に向上した。さらに、自由表面に厳密な応力境界条件を課すことができ、自由表面形状を正確に予測できる

2.2

時間進行法と空間離散化

移動サーフェスメッシュ LSMPS 法 へと発展しており、時間進行法は松永ら（6）によって提案された新たなフラ

この移動サーフェスメッシュ LSMPS 法を用いることで、クショナルステップ法を用いる。この手法では、（2）式

未だ確立されていない粘弾性流体の解析を高精度に再の対流項を陽的オイラー法にて定式化（（6）式）し、粒子現・予測することが期待できる。

移動による速度と応力変化を Taylor 展開近似によって（4）

そこで本研究では、移動サーフェスメッシュ LSMPS法に粘弾性モデルを適用し、粘弾性流体の挙動を高精度求める（（11）、（12）式）。
具体的な計算手順として、まずは仮速度 * を求める。に予測する演算手法を新たに開発したので、その内容に（6）ついて述べる。

2.

演算アルゴリズム

添字 は時間ステップ、Δt は刻み時間である。圧力は圧力ポアソン方程式を陰的に解いて求める。

2.1

支配方程式

（7）

非圧縮性流体の支配方程式は連続の式とナビエストークス方程式、および粘弾性流体を表す構成方程式としてOldroyd-B モデル（5）を用いる。

次に、現在位置（1）テンソルでの新しい時刻の速度と応力を計算する。
（8）

（2）
（9）

（3）
最後に粒子を移動し、新しい位置は流速ベクトル、ρは密度、 は時間、 は圧力、μs は応力テンソルでの速度とを求める。

溶媒粘性率、τは応力テンソル、 は重力ベクトル、λは緩和時間、μ は高分子粘性率、（10）は変形速度テンソルである。本手法の特徴として、（2）式に対流項を含む点にある。ラグランジュ法では粒子（計算点）が移動す（11）

ることで対流項の役割を果たすため、対流項を計算式に含める必要がないと考えられており、従来の MPS 法とSPH 法も対流項は現れない。この対流項の取扱いが、タ（12）イムステップごとに圧力誤差が発生する要因とされており、本手法では対流項を考慮することで、誤差の発生を回避でき計算精度の向上が見込める。

は近似の次数を決める正の整数、自由表面では次式に示す応力の釣り合いを境界条件とは位置の修正距離である。位置の修正は粒子配置が乱雑になっていくのを防ぐために必要な処理となる。タイムステップごとにして与える。

流体粒子が移動する際、計算時間の経過とともに粒子の（4）は流体周囲の気体圧力、 は外向き単位法線ベクト
ル、（2）は表面張力ベクトル、は気体粘度である。

乱雑さが増していき、最終的に数値不安定性や計算誤差,。に繋がるからである（7）（8）

流体内部粒子の空間離散化には、重み付き最小二乗法を用いた LSMPS スキーム（3）を用いる。標準 MPS 法は日本製鋼所技報 No.74（2023.11）

技 術 論 文

MPS 法による粘弾性流体の数値シミュレーション

勾配（∇）とラプラシアン（∇2）の微分演算子に対して、あは未知数ベクトル、は重み関数、は影響半径でる粒子 i とその近傍の粒子 j との間で粒子間相互作用モある。重み付き二乗誤差 J を最小化する正規方程式は次デルを与え離散化を行うが、この粒子間相互作用モデルのようになる。

は等方的な粒子配置を仮定し導出しているため、計算過程中に粒子が不規則な配置になると予測精度が著し（22）く低くなることが知られている。一方で LSMPS 法はこの問題を改善するために考案された手法で、Taylor 展（23）

開と最小二乗法を用いて空間離散化を施すことで、不規則な粒子配置に関わらず高精度の予測が可能となる。

が特異行列でない場合はただ一つのLSMPS 法には type-A と type-B の形式があるが、本論

モーメント行列文では、最も基本的な LSMPS type-A を用いた離散化解が存在し、正規方程式（22）式の解が得られる。

について説明する。
ある粒子 i の近傍粒子 j における関数φの Taylor 展開（24）は以下となる。
（13）

最終的に、勾配（∇）とラプラシアン（∇2 ）は以下のように計算することができる。

（14）
（25）
であり、は近似の次数を決める正の整数、（26）はスケーリング係数である。（14）式は以下のように表すことができる。

3.

（15）
は多項式ベクトル、はスケーリング行列、は微分作用素ベクトルである。例えば、 = 2 の場合は以

3.1

検証解析

円筒回転問題による検証解析

本解析手法の妥当性を検証するために、図 1 に示す円筒回転問題での解析を行った。円筒半径 r0 = 0.1 m の内下のように成分表示される。

部に高さ h0 = 0.2 m の流体を配置し、重力加速度（16）= 10m/s を与え回転速度Ω = 5, 10 rad/s の 2 条件で回転さ2せた。流体は密度ρ= 1000 kg/m3、粘度η= 0.001 Pa・s
（17）（18）の純粘性ニュートン流体とし、粒子間距離（粒子径）l0を 1.25 × 10 -3 m と 2.5 × 10 -3 m の 2 流体とした。また、解析対象領域は 2 次元軸対称問題として扱った。

x ij = x j - x i 、y ij = y j-y i である。（15）式において、重み付き二乗誤差 J は以下となる。
（19）
（20）

（21）
図 1 円筒回転問題の計算モデル（3）

技 術 論 文

MPS 法による粘弾性流体の数値シミュレーション

図 2 に粒子間距離 l0 = 1.25 × 10-3 m、回転速度Ω =10 rad/s とした条件での時間ごとの圧力分布と速度分布の予測結果を示す。圧力分布の結果から、流体は重力の作用によって自由表面から底面に向かって圧力が高くなっており、妥当な予測結果が得られている。また、どの時間に対しても圧力振動が発生しておらず、滑らかな圧力分布が得られている。これは粒子配置の乱れに応じた空間離散化が行えているためであり、LSMPS の予測精度の高さが窺える。速度分布については、計算開始から 3秒にかけて、流体粒子が外側円筒に向かうような速度ベクトルとなっており、これは外側円筒の回転によって流

図 3 自由表面形状の計算結果と理論解との比較体に遠心力が作用しているためと考えられる。その後 4秒から次第に速度が小さくなり、5 秒では半径と高さ方

3.2

ワイゼンベルグ効果シミュレーション向の速度はほぼ 0 と平衡状態に達していることが分かる。

同心二重円筒間に一定体積の流体を配置し、内側円筒次に、自由表面形状の理論解と予測結果との比較評価（以下、中心ロッド）を回転させワイゼンベルグ効果のシミュレーションを行った。座標系は 2 次元軸対称とし、を行った。次式に理論解を示す。

シミュレーションモデルを図 4 に示す。計算条件は、粒（27）子間距離 l 0 = 1.25 × 10 -3 m、溶媒粘度比率β =（μs ⁄（μs + μ ））= 0.5 、回転速度Ω = 5 rad/s とし、緩和時は半径位置 r での界面高さ、 は初期界面高さである。図 3 に各解析条件での予測結果とそれに伴った理間λ = 0.1, 0.3, 0.5 s の 3 ケースでの粘弾性流体、およと回転速度度場の計算結果を示す。λ = 0.1 s の場合では、中心ロッΩを変化させた全 4 ケースとも予測結果と理論解とでド近傍で下向きの速度場になっており、液面が低下し良好な一致が見られ、本手法は高精度の予測が可能でている。一方λ = 0.5 s では、中心ロッド近傍で上向きの速ある。

度が発生し液面 が 盛 り 上がっていることが確認でき る。

論解を示す。自由表面形状は、粒子間距離びニュートン流体とした。図 5 にλ = 0.1, 0.5 s での速図 2 円筒回転問題の計算結果（4）

日本製鋼所技報 No.74（2023.11）

技 術 論 文

MPS 法による粘弾性流体の数値シミュレーション

図 6 に各λの条件とニュートン流体の定常状態での自由側方向への力が大きくなる。その結果、外側の流体は内表面形状の計算結果、およびニュートン流体の理論解を側へと流れ、中心ロッド部で盛り上がる挙動となる。

示す。ニュートン流体は外向きに遠心力が働くことで、次に、粒子間距離 l0 = 6.25 × 10 -4 m の高解像度とし、中心ロッドから流体が離れる挙動となり、中心ロッド近回転速度Ωを変化させて計算を行った際の、流体がロッ傍で液面が初期液面高さよりも低くなっている。また、ドに巻き付く高さ h について理論解（9）と比較した。以下

自由表面形状は理論解と良好に一致しており、計算結果に理論解を示す。の妥当性が示唆された。粘弾性流体では、λが大きくなるにつれ、中心ロッド近傍での液面が高くなる結果が得（28）られた。中心ロッドが回転すると、内側の流体に生じる応力は外側に比べ大きく、また伸びも大きくなる。粘弾性流体は伸びが大きいほど、元の状態を維持する張力、（29）すなわち法線応力差が発生し、この法線応力差による内（30）R は中心ロッド半径、σは表面張力、ρは流体密度、Wi はワイゼンベルグ数、ξはパラメータである。図 7に各回転速度Ωでの h ⁄ R と Wi 2 の関係について、予測結果と理論解をそれぞれ示す。予測結果は、ワイゼンベルグ数が小さい（回転速度が小さい）場合に理論解とよく一致しており、一方で、回転速度が大きくなるにつれて予測結果が理論解から徐々に乖離している。この原因
は、巻付き高さの理論解（28）式の右辺第 2 項において 4
図 4 ワイゼンベルグ効果の計算モデル次以上の項を無視しており、回転速度が増すにつれて理論解自体の誤差が次第に大きくなるためと考えられ、これらの誤差を考慮すると、計算結果は理論解と同等の結果が得られている。
以上のことから、本解析手法によって粘弾性流体のワイゼンベルグ効果を高精度に再現できている。

（a）λ = 0.1 s（b）λ = 0.5 s

図 5 速度分布の計算結果

図 7 巻付き高さの計算結果と理論解との比較

図 6 粘弾性流体とニュートン流体の自由表面形状、および理論解（5）

技 術 論 文

MPS 法による粘弾性流体の数値シミュレーション

4.

おわりにmethod”, International Journal for Numerical
Methods in Engineering, Vol. 121（2020）, pp. 4101-

移動サーフェスメッシュ LSMPS 法に粘弾性流体の4133

Oldroyd-B モデルを適用し、粘弾性流体の挙動を高精度（9） G. S. Beavers and D. D. Joseph:“The rotating rodに予測する演算手法を新たに開発した。本解析手法にviscometer”, Journal of Fluid Mechanics, Vol. 69

よって粘弾性流体で発生するワイゼンベルグ効果を再（1975）  No. 3, pp. 475-511

現・予測でき、自由表面形状の理論解ともよく一致した。
以上より、本報の解析手法は粘弾性流体の挙動を予測できる有用な手法といえる。今後は、ダイスウェル現象などの粘弾性流体で生じる特異的な現象の予測や、非線形粘弾性モデルである Giesekus モデルや PTT モデルへの拡張を行う所存である。

参考文

献（1） S . Koshizuka and Y. Oka : “Moving particles em i - i mpl icit met hod for f rag ment at ion ofincompressible f luid ”, Nucl. Sci. Eng., Vol. 123（1996）, pp. 421-434（2） J. Monaghan:“Smoothed particle hydrodynamics”,Ann. Rev. Astron. Astrophys., Vol. 30（1992）, pp.
543-574（3） T. Tamai and S . Koshizuka : “Least squares
m o v i n g p a r t i c l e s e m i - i m p l i c i t m e t h o d ”,
Computational Particle Mechanics, Vol. 1（2014）  
No. 3, pp. 277-305（4） T. Matsunaga, S. Koshizuka, T. Hosaka and E.
Ishii:“Moving surface mesh-incorporated particle
method for numerical simulation of a liquiddroplet”, Journal of Computational Physics, Vol.
409（2020）, p. 109349（5） J. G. Oldroyd:“On the formulation of rheologicalequations of state”, Proceedings of the Royal
Society of London. Series A. Mathematical andPhysical Sciences, Vol. 200（1950）No. 1063, pp.
523-541（6） 松永拓也 , 越塚誠一 :“粒子法における時間進行法の改良”, 日本機械学会論文集 , Vol. 87（2021）  No.
895, 20–00437（7） R. Xu, P. Stansby and D. Laurence:“Accuracyand stability in incompressible SPH （ISPH）
based on t he projection met hod and a newapproach”, Journal of computational Physics, Vol.
228（2009）   No. 18, pp. 6703-6725（8） T. Matsunaga, N. Yuhashi, K. Shibata and S.
Koshizuka:“A wall boundary treatment using
analytical volume integrations in a particle（6）

日本製鋼所技報 No.74（2023.11）