弾性流体潤滑におけるキャビテーションに関する研 究

九州大学学術情報リポジトリ

Kyushu University Institutional Repository

弾性流体潤滑におけるキャビテーションに関する研 究

大津, 健史

九州大学大学院工学府

https://doi.org/10.15017/21997

出版情報：Kyushu University, 2011, 博士（工学）, 課程博士 バージョン：

権利関係：

第 3 章 キャビティー成長のモデル

3.1 Initial stageのキャビティー成長モデル

本節では，第2章の実験結果において分かったキャビティー成長現象について考察し，成長のモデ ルを検討する．ここでは，特に，第2章2.2に示した雰囲気気体の影響に関する実験結果を基に考察，

検討を行う．

2.2 に示した実験結果より，キャビティー長さの時間依存性と雰囲気依存性が確認された．また，

発生からの成長過程は2段階(Initial stage，Second stage)であることも確認された．発生後すぐに起こ る成長過程である Initial stage では，キャビティーは急速に成長し，溶解気体の影響が見られない．

Initial stage後に起こる成長過程であるSecond stageでは，キャビティーの成長速度は低下し，溶解気

体の影響が見られる．また，キャビティーは，溶解度の高い気体ほど長くなる．これらの結果を基に，

本節ではInitial stageのキャビティー成長のモデルを検討する．

3.1.1 キャビティーモデルの提案

(1) キャビティー成長速度とキャビティー形状からの考察

Fig. 3.1-1に，発生から11 msにおけるPAO63(動粘度63mm²/s@313K)でのキャビティー，8 msにお

けるPAO400(動粘度397mm²/s@313K)でのキャビティーをそれぞれ示す．各滑り速度は，11，20 mm/s

であった．この図より，Initial stageにおけるキャビティー形状は，Fig.2.2-10にも示したように，急 速な成長を示した特異的なものであることが確認できる．キャビティー後端は，一部分が伸び，複数 の突起をもった形状となっている．この形状は， Second stageにおけるキャビティーのように滑らか な界面を持ったものとは異なる．

Dellis, P. and Arcoumanis, C. ¹は，ピストンリング－シリンダライナに発生するキャビテーションを

観察し，その成長について報告している．Fig.3.1-2に，クランク角20°，90°でのキャビティーの画 像と油膜圧力の測定結果を示す．この図において，20°のときに発生しているキャビティーは，全体 が一様に拡大するのではなく，部分的に伸び，シダの葉のような形状を示している．この形状は，

Fig.3.1-1に示すキャビティー形状と似ている．また，この時の油膜圧力は，絶対真空圧まで低下して

いる．従って，この結果より，キャビティーの発生とその直後の成長は，絶対真空まで低下した負圧 の発生と起因していると考えられる．その後，90°において，キャビティーは成長し，潤滑油のメニ スカス領域を越え，大気とつながっている．また，潤滑油が筋状に流れていることも確認できる．油 膜内圧力は，大気圧と同等である．

また，部分的に成長する界面を示す相変化の現象として，デンドライト結晶成長が知られている^2,3．

Fig.3.1-3に，デンドライト結晶成長の画像を示す．図に示すように，一部分で伸びた，複数個の突起

が界面には見られる．Fig.3.1-4にて，その成長機構を説明する．融液が凝固していく過程((a)t = t1)に おいて，液相側が凝固点よりも低い温度であると(過冷却状態)，一部分の液体が凝固し突起のように なる(図中のA部分)．さらに，時間t = t2において，突起の前方の位置(A)でも，凝固点よりも低い温 度の過冷却状態が起こっているならば，凝固が継続され，突起が成長していく((b) t=t2)．突起の両側 の液相には，固相となった突起領域からの溶質原子の排出があり，凝固点が低下する(凝固点降下現 象)．そのため，両側の部分は，液体状態を維持する．従って，時間経過とともに，突起部分(A部分) のみが凝固し，伸びていき，図に示すような形状となる((c) t=t3)．また，Fig.3.1-5(四臭化炭素，CBr₄ における結晶成長)に示されているよう^3,4に，温度勾配がx方向に凝固点温度以下の温度を維持する 条件(図中の条件(a))ならば，結晶成長は速く，突起の数も多くなる．一方で温度勾配が凝固点温度よ りも大きな条件(図中の条件(b))であるならば，成長速度は遅く，突起の数も少なくなる．

これらのことを考慮し，Initial stageにおけるキャビティーの成長機構を考えると以下のようになる．

開始後，接触域後方では，溶解気体の析出圧力を越える負圧が発生し，キャビティーが形成される．

その後も，キャビティー後方には大きな負圧が発生していれば，直前まで大気圧条件下にあった溶解 気体は，次の瞬間に過飽和の状態となり，液体中には溶解できなくなる．その結果，溶解気体は，キ ャビティー内への析出する．時間とともに，この析出が起こるため，キャビティーが成長すると考え られる．

(2) 気泡力学からの考察

気泡力学では，気泡径の変化が，気泡内の圧力と液体の圧力差によって起こるとされている．以下 の式は，Rayleighの式⁵で，気泡径Rの気泡の運動方程式を示している．

) 1 (

2 3 ₂

pb

p R

R

R + = − _∞+

&

ρ

&

& (3.1.1-1)

R：気泡径，ρ^{：密度，p}∞：液体中の無限遠方での圧力，p_b：気泡内の圧力

また，表面張力を考慮した，Rayleigh-Plessetの式は，以下のように示される．



 



− + −

=

+ _∞

p R p R

R

R _b

σ

ρ

2 1

2

3 &₂

&

& (3.1.1-2)

σ: 気液間の表面張力

なお，表面張力と気泡径の関係は，次の Young-Laplace の式で考えられている．この式において，

気泡内圧力は，液体の圧力よりも表面張力の分だけ高くなる．

p R p_{b w} 2

σ

=

− ^(3.1.1-3)

pw: 液体の圧力

気泡力学の視点から考えると，以下のようなキャビティーの成長機構が考えられる．

潤滑油中には，溶解気体分子が凝集した気泡核が存在するといわれている^6,7．大気圧下では，この 気泡核内の圧力は，周囲と同じ大気圧，あるいは，式(3.1.1-3)で示されるように大気圧よりも高くな る．気泡核が潤滑油とともに移動し，接触域後方の負圧領域に進入すると，気泡内と周囲の圧力差に よって，式(3.1.1-1)で示されるように，気泡核は自身の領域(径)を拡大する．その圧力差が大きいほど，

この拡大は速くなり，拡大する領域は大きくなると考えられる．従って，キャビティーは急速に成長 する．

（1），(2)において示した両モデルともに，接触域後方の圧力分布を調べ，それを基に検討する必 要がある．次に，圧力分布を調べた結果を示す．

3.1.2 接触域後方の圧力分布の計算

（1） 滑り試験における圧力分布の計算

3.1.1の(1)，(2)のモデルの提案で述べたように，Initial stageにおけるキャビティーの発生には，接

触域後方に発生する負圧分布が関係していると考えられる．ここでは，その圧力分布を，レイノルズ 方程式を解くことにより計算し，キャビティーの成長との関係を検討する．

Fig.3.1-6 に，本計算のモデル図を示す．本モデルは接触域の後方部分におけるすきまを取り扱い，

その領域において, 以下の 1次元 x方向のレイノルズ方程式(式(3.1.2-1))を解いた．式(3.1.2-2)は，式

(3.1.2-1)をxで積分したものである．計算範囲は，Fig.3.1-6に示されるように，x1からメニスカスの半

径位置x2までとした．x1で定義される位置は，ヘルツ接触半径aとキャビティー長さxcの和で示され る．メニスカスの半径位置x2は，実験において観察された5 mmと設定した．また，本計算では，境 界条件として，キャビティーと潤滑油の境界での圧力P_x1＝0(大気圧)，メニスカス位置での圧力P_x2

＝0(大気圧)と仮定した．

dx dh U dx dp h dx

d

2 12

3 =





 η 6 3

h h U h dx

dp =

η

− _m ^(3.1.2-2)

p：圧力，η：潤滑油粘度，U：滑り速度，h：すきま，hm：積分定数

すきま形状h(x)は，ヘルツ接触を考慮した以下の式⁸を用いて計算した．

 





 





 

 



 −

 +



 



 −

− +

=

⁻

5 . 0

2 2 1

2 2 2

min

2 cos 1

)

( a

x x a a

x R

h a x

h π

^(3.1.2-3)

a：ヘルツ接触円半径，R：レンズの曲率半径，h_min：最小油膜厚さ

また，h_mは，以下のように示される．

∫

=

∫

2 1

2 1

3 2

1 1

x x x x m

h dx h dx h

圧力分布は，式(3.1.2-2)を積分し，以下のように示される．



 



 −

=

∫

m

∫

x^x

x

x dx

h h h dx U x

p

1

1 2 3

1 6 1

)

(

η

(3.1.2-1)

(3.1.2-4)

(3.1.2-5)

Fig.3.1-7に，計算された圧力分布の時間変化を示す．計算に使用する各時間におけるキャビティー 長さは，PAO63の空気中での結果 (Fig.2.2-11) を使用した．また，計算に使用した滑り速度，潤滑油 粘度等は，Fig.2.2-11 の実験のものと同じ値である．なお，この図において，圧力はゲージ圧で示し ている．

この図より，0 msに，絶対真空の負圧の発生が見られ，3.4 msまで絶対真空圧の発生が続く．その 後，時間とともに，発生する負圧の大きさは小さくなり，緩和していくことが分かる．Fig.3.1-8 に，

キャビティー長さの変化とFig.3.1-7に示される圧力の最小値の時間変化を示す．この図より，3.4 ms までの時間では，絶対真空圧の発生とキャビティーの急速な成長が関係していることが分かる．その 後，時間とともに発生する圧力が高まる，このとき，キャビティーの成長速度は低下していることも 分かる．

従って，キャビティー後方に発生する負圧とその緩和が，Initial stageにおけるキャビティー成長に 影響を与えていると考えられる．

(2) 引き離し試験における圧力分布の計算

Fig.3.1-9に，計算モデルを示す．計算は，レンズが完全にディスクから離れた時間である 2 msか

ら行った．なお，r1は，キャビティー径の半分の値の位置である．

レイノルズ方程式は，以下の式(3.1.2-6)である．なお，接触域は円形であるので，レイノルズ方程 式は，極座標系のもの⁹を使用した．式(3.1.2-7)は，式(3.1.2-6)をr で積分したものである．境界条件 として，r₁の位置での圧力をP_r1=0(大気圧)，メニスカス位置r2での圧力をP_r2=0(大気圧)とした．

r2は，滑り試験時の計算と同じく，5 mmとした．

dt v dh dr h dp dr

d = =



 



 3

12 1

η

(3.1.2-6)

1 3 3

12 12

C h h vr

dr

dp η η

+

= (3.1.2-7)

v ( = dh/dt)：引き離し速度，C₁：積分定数

すきま形状h(x)は，円の方程式より，以下の式で計算される．

(

^{2 2}

)

) 0

(r h R R r

h = + − − (3.1.2-8) h0: 中心位置における膜厚，R: レンズの曲率半径

式(3.1.2-7)を，境界条件を使用し解くと，以下のようになる．

∫

∫

⁺

= ^r

r r

r dr

C h h vrdr

p

1

1 3 1 3

12

12η η

(3.1.2-9)

ここで，C₁は，以下のように示される．

∫

−

∫

=

2 1 2 1

3 3

1 1

1

r r r r

h dr h vrdr

C (3.1.2-10)

計算に使用するキャビティー径は，Fig.2.2-21 における空気中での実験結果を使用した．また，引 き離し速度，潤滑油粘度は，実験のものと同じ値である．

また，引き離し過程では，式(3.1.2-8)中の膜厚 h0が時間とともに変化する．従って，計算時間にお ける各時間で，膜厚h0を設定することになる．本計算では，2，3 msにおいては，Fig.2.2-19から読 み取り，それ以降の時間においては，引き離し速度によって求めた値を使用した．

Fig.3.1-10 に，計算された圧力分布の時間変化を示す．なお，この図において，圧力はゲージ圧で

示している．

2，3 msにおいて，絶対真空の圧力が発生しており，この負圧の発生が Fig.2.2-20におけるキャビ

ティーの発生と急速な成長と関係していると考えられる．その後，発生する圧力は上昇している．一 方で，実験結果では，この時間においても，キャビティーの成長速度は大きい(Fig.2.2-21)．これは，

計算条件の膜厚h0の設定による影響と考えられる．時間4 ms以降は，膜厚の設定を引き離し速度か ら求めているが，実際の膜厚は，この計算値よりも小さいことが考えられる．

3.1.3 キャビティー成長モデル

3.1.2 における負圧分布の計算により，Initial stageにおけるキャビティー成長は，負圧の発生と関

係していることが示唆された．ここでは，3.1.1 で提案したキャビティーの成長モデルを再度検討す る．

まず，3.1.1(1)で提案したモデルについて検討する．

Fig.3.1-11には，時間0 sにおける，接触域後方に発生する負圧分布を示す．この図では，絶対真空

圧の処理を行っていない．接触域出口から約3 µm後方の位置で，最小圧力‐55MPaの負圧が発生し ている．液体は，負圧が加わったときに，その張力に耐え，液体を維持することができると報告され ている¹⁰．パラフィン系鉱油では，約12 MPaの張力に耐えることができると報告されている¹¹．よ って，負圧が発生している系においても，液体はその負圧によって加わる張力に耐え，破断しないこ とが考えられる．本計算で求められた負圧は，-55MPaで，その張力(例えば，12 MPa)と比較して高い．

従って，潤滑油分子は破断し，空洞が形成されると予想できる．これがキャビテーション核となる．

なお，潤滑油中に気泡核などが存在すると，キャビテーション発生の圧力は上昇する ^7,12．また，せ ん断場では，液体の張力が低下するため，キャビテーション核の形成が容易になると考えられている (これは”Triboneucleation”と呼ばれる)¹³．

0 sの時点で，負圧の発生によりキャビテーション核が形成される(Fig.3.1-12)．ディスクが移動する とき，このキャビティー核もディスクとともに移動すると考える．接触域出口からは，油膜厚さ分の 油量のみ供給されるので，キャビティー領域は拡大していくことになる．ただし，接触域後方では，

負圧が発生しているために，圧力流れにより，潤滑油は接触域の方へ逆流する．従って，この逆流に よって，キャビティー領域は消滅することが考えられる．しかし，実験結果では，キャビティー領域 は消滅することなく拡大を続け，発生初期のキャビティーの成長速度は滑り速度に近づいている．従 って，油の逆流の影響を受けずに，キャビティー領域はディスクの移動とともに拡大していることが 示唆される．これは，次のような溶解気体の析出が関係していると推測される．

Fig.3.1-13には，0 s直後の時間における接触域後方を模式的に示した図である．拡大したキャビテ

ィー領域の後方には，絶対真空圧を越える負圧が発生している．その部分は，直前まで大気圧であっ たが，この時間では圧力が絶対真空圧まで低下している．潤滑油中の溶解気体濃度Cは，以下のヘン リーの式(3.1.3-1)で示されるように，気体の分圧(圧力)Pと比例する．大気圧P0での濃度C0＝P₀/Hは，

負圧Pnが発生した領域での濃度Cn＝Pn/Hよりも高い(C0>Cn)．従って，負圧が発生した直後に，溶解 気体濃度は，C₀から Cnに変化しなければならなくなる．このとき，溶解気体は液体中に過飽和状態 となり，エネルギーが高くなる．

HC

P =

(3.1.3-1)

P：圧力, Pa，H：ヘンリー定数, Pa/(mol/m³)，C：溶解濃度, mol/m³

過飽和となった溶解気体は，キャビティー領域と潤滑油の界面から，その分圧差(濃度差)によって，

キャビティー側へ析出していく．析出することにより，キャビティー内には気体が存在するので，キ ャビティー内の圧力が上昇し，周囲の潤滑油の圧力と圧力バランスをとる．従って，このキャビティ ーは，安定的に潤滑油中で存在できるようになる．このように，キャビティー領域には潤滑油が流入 するより先に溶解気体が析出したため，その領域を拡大できたと考えられる．

Fig.3.1-14には，Fig.2.2-10に示す10 msのときの接触域後方の模式図である．この時間では，キャ

ビティーの急速な成長は起こっていない．接触域後方に発生する負圧は，Fig.3.1-13 のときと比較し て，緩和している．従って，キャビティー領域後方の潤滑油中の溶解気体量の過飽和度は小さくなっ ている．そのため，キャビティー内への溶解気体の析出量は小さくなる．また，この時間では，潤滑 油の逆流が，キャビティー領域の拡大に対して抵抗となるため，急速な拡大はなくなる．キャビティ ー領域は，気体の析出により内部の圧力が上昇しているため，この時間においても，安定的に潤滑油 中に存在できる．

従って，この時間域でのキャビティーの成長は，以下のようにまとめられる．

負圧によって，接触域出口でキャビティー核が形成され，ディスクの移動とともに，その領域の拡 大が起こる．それとともに，キャビティー領域後方では，負圧の発生により過飽和となった溶解気体 のキャビティー内への析出が起こる．そのため，キャビティー内圧力が上昇し，潤滑油中に安定的に 存在できるようになる．時間とともにキャビティー領域後方に発生する負圧は緩和し，溶解気体の析 出量が小さくなる，また，後方での潤滑油の逆流も影響するために，次第にキャビティーの成長速度 が小さくなる．

実験結果より，キャビティーの成長に，雰囲気の影響は見られない．しかし，キャビティー内圧力 においては，雰囲気の差がみられると考えられる．Fig.2.2-21 より，キャビティーが収縮する過程に おいて，溶解気体の影響が見られ，収縮後の気泡径は，ヘリウム，アルゴン，二酸化炭素の順に小さ い．収縮過程におけるキャビティー部への潤滑油の流入は，潤滑油とキャビティー内の圧力差による 圧力流れに起因したものである．従って，二酸化炭素中では，キャビティー内圧力が高いため，キャ ビティー領域への流入量が小さくなり，収縮に時間がかかると考えられ，ヘリウム中では，キャビテ ィー内圧力が低いため，流入量が多くなり，収縮が早いと考えられる．この雰囲気によるキャビティ ー内圧力の違いは，溶解気体量に依存した結果であると考えられる．従って，二酸化炭素中では，溶 解気体量が多いため，キャビティー領域に析出する気体量が多くなり，キャビティー内圧力が上昇す ると推測される．一方で，ヘリウム中では，溶解気体量が少ないため，キャビティー領域へ析出する 量が少なく，キャビティー内圧力が低いと推測される．キャビティー領域の拡大は負圧の影響を受け

るため，雰囲気の大きな影響は見られないが，キャビティー内に析出する気体の量，および，キャビ ティー内の圧力は雰囲気によって異なることが分かる．

次に，3.1.1(2)のモデルについて検討する．

このモデルでは，負圧部に進入した気泡核の径の拡大によって起こるものである．従って，Fig.3.1-15 のような機構が考えられる．

潤滑油中の気泡核は，負圧部に進入することにより，式(3.1.1-1)で示されるように，圧力変化によ って領域を拡大する．図のように，キャビティー後方に大きな負圧の発生する時間，図におけるt=t1

では，その負圧部分に進入した気泡核が拡大し，体積が V1から V2へ変化する．その後，t=t₂でも，

負圧の発生により，気泡核は，V₃から V4 へと体積を変化させる．それぞれが合体することにより領 域が拡大していく．時間経過とともに負圧の発生が緩和すると，圧力変化が小さくなるので，気泡核 の拡大も小さくなる．従って，キャビティーの成長は小さくなる．

このモデルでは，溶解気体量と気泡核数の関係が強く影響することが考えられる．従って，溶解気 体量が多いと，気泡核も多くなり，キャビティーの成長速度や形状が変化することも予想される．気 泡核の測定などの測定は，現在までに行っておらず，この点における影響は不明である．この点に関 する検討は今後の課題としたい．

3.1.4 まとめ

本節にて，Initial stageの各成長過程において，その成長のモデルを検討した．

キャビティーの急速な成長は，接触域後方に発生する負圧領域と関係がある．また，成長速度の低 下は，負圧の緩和と関係がある．従って，時間とともに変化する負圧領域の変化を把握することで，

この段階におけるキャビティーの成長を理解できる．

この結果を基に，以下のモデルが提案される．

負圧の発生とともにキャビティー核が形成され，その拡大が起こる．このとき，過飽和となった溶 解気体がキャビティー内へ析出するために，キャビティーが安定的に存在できるようになる．負圧が 緩和すると，溶解気体の急速な析出がなくなり，キャビティーの成長が小さくなる．

また，潤滑油中の気泡核が負圧領域に進入するとともに拡大し，キャビティー領域を拡大させる．

負圧が緩和すると，気泡核の拡大が小さくなり，成長も小さくなる．

(参考文献)

1 Dellis, P. and Arcoumanis, C., “Cavitation development in the lubricant film of a reciprocating piston-ring assembly,” Proceedings of Institution of Mechanical Engineers, Part J, Journal of Engineering Tribology, 218, 2004, 157-171

2 Kurz, W. and Fisher, D.J., Fundamentals of Solidification, Trans Tech Publications, 1998 3 後藤, 結晶成長, 内田老鶴圃, 2003

4 Eshelman, M.A., Seetharaman, U.S. and Trivedi, R., “CELLULAR SPACING –Ⅰ. STEADY STATE GROWTH”, Acta Metllurgica, 36, 4, 1988, 1165-1174

5 日本流体力学会, 混相流体の力学, 朝倉書店, 1991

6 中原, “流体油膜の挙動 -キャビテーション-“, 潤滑, 26, 3, 1981, 146-152

7 松本, 奥平, 和田, 榎本, 市川, ”気泡核分布に及ぼすキャビテーションの影響”, 日本機械学会論文 集(B編), 51, 472, 3844-3851

8 Wedeven, L.D., Evens, D. and Cameron, A., “Optical Analysis of Ball Bearing Starvation,”

Transaction of the ASME, Journal of Lubrication Technology, 93, 1971, 349-363 9 Cameron, A., Principles of Lubrication, John Wiley and Sons INC, 1966

10 中井, 沖野, “潤滑油膜中の張力”, 日本機械学会論文集(第3部), 40, 333, 1974, 1465-1471

11 Vincent, R.S. and Simmonds, G.H., “Examination of the berthelot method of measuring tension in liquids”, Proceedings of the Physical Society, 55, 1943, 41-48

12 Knapp, R.T., “Cavitation and nuclei”, Transaction of ASME, 80, 6, 1958, 1315-1324

13 Hayward, A.T.J., “Tribonucleation of bubbles”, British Journal of Applied Physics, 18, 5, 1967, 641

(a) PAO63, 11 mm/s (Cavity at 11 ms)

100 µm

(Cavity at 8 ms) (b) PAO400, 20 mm/s

Fig.3.1-1 Cavity shape in growth of initial stage

Fig.3.1-2 Relation between cavity shape and oil film pressure ¹

Solid

Liquid

(a) t = t₁

Fig.3.1-4 Schematic figure of dendrite crystal growth ³ (b) t = t₂ (c) t = t₃ Fig.3.1-3 Dendrite crystal growth ²

Solid

T

T Ts

Ts

Ts: Freezing temperature

Solid Liquid Liquid (b)

(a)

Fig.3.1-5 Effect of temperature gradient on growth rate (in CBr4) ³

x

x

Fig.3.1-6 Calculation model in sliding test

Lens

Disc

h(x)

U

x x ₁ = a + x

_c

x 2

h

min

Cavity length x

_C

a x ₁

Fig.3.1-7 Changes in pressure distribution with time

0 100 200 300 400 500

-1 -0.8 -0.6 -0.4 -0.2 0

P re s s u re , 1 0 5 P a

Position, µm 0ms

12.3ms 14.6ms

1.1ms 2.2ms 3.4ms 4.5ms 5.6ms 6.7ms 7.8ms 10.1ms

0 5 10 15

-1 -0.8 -0.6 -0.4 -0.2 0

0 50 100 150 200

Time, ms

P re s s u re , 1 0

5

P a C a v it y l e n g th , µ m

Pressure Cavity length

Fig.3.1-8 Relation between cavity growth and minimum pressure

0 500 1000 -1

-0.8 -0.6 -0.4 -0.2 0

P re s s u re , 1 0

5

P a

Position, µm

2ms 3ms 4ms 5ms 6ms

Fig.3.1-9 Calculation model in separation test

Fig.3.1-10 Changes in pressure distribution with

r

r ₂ h(r)

h ₀

v

O R

Cavity size 2r ₁ Lens

Disc

r

r ₂ h(r)

h ₀

v

O R

Cavity size 2r ₁ Lens

Disc

Fig.3.1-11 Pressure distribution at 0 s -60

-50 -40 -30 -20 -10 0

0 100 200 300 400 500

Position, µm Pressure, ×106 Pa

Oil Disc

Lens

P re s s u re

-0.1MPa

0

Oil Disc

Lens

P re s s u re

-0.1MPa

0

Oil Disc

Lens

P re s s u re

-0.1MPa 0

Oil Disc

Lens

P re s s u re

-0.1MPa 0

Dissolved gas

Cavity

Dissolved gas

Cavity

Fig.3.1-12 Initiation of cavity at 0 s

P re ssu re P re ssu re

Oil Disc Lens

P re s s u re

-0.1MPa 0 Bubble core

t＝t₁

V

1

t＝t₂

V

2

V

3

V

4

Oil Disc Lens

P re s s u re

-0.1MPa 0 Bubble core

t＝t₁

V

1

t＝t₂

V

2

V

3

V

4

Fig.3.1-14 Cavity growth after initial stage

Oil Disc

Lens

P re s s u re

-0.1MPa 0

Oil Disc

Lens

P re s s u re

-0.1MPa 0

Dissolved gas

Cavity

Cavity

P re ssu re

3.2 MDシミュレーションによるキャビティー発生の検討

本節では，分子動力学法を使用した数値計算を行い，負圧が発生した液体中でのキャビティーの発 生，および，それに対する溶解気体の影響を調べた．

3.1で，Initial stageにおける成長モデルは，キャビティー後方に発生している負圧領域と関係して いることを示した．その一つのモデルとして，負圧の発生によるキャビティー核の発生と，その拡大 とともに溶解気体がキャビティー内へ析出するために，キャビティーが成長していくと述べた．この モデルの中で，負圧が発生した直後に起こるキャビティーの形成と，キャビティー発生後に起こる溶 解気体の析出現象は，その規模と時間スケールから，実験では確認できない．従って，本節では，分 子動力学法を用いたシミュレーションを行い，その現象を調べ，その結果より提案したモデルを検討 した．これにより，Initial stageにおけるキャビティー成長がより理解できる．

本計算では，液体に水，溶解気体にヘリウム，アルゴンを使用した．上下の鉄壁面に水と溶解した 気体を挟み，上部壁面を移動させることにより，系の圧力を低下させた．この結果より，圧力低下に より発生するキャビティー，および，このキャビティーの発生に対する溶解気体の影響を検討した．

3.2.1 分子動力学法

分子動力学法(Molecular Dynamics Simulation, MD法)は，物質を分子の集合体として取り扱い，一つ 一つの分子の動きを運動方程式から求め，その時間における分子の一連の動きを追求し，それらを基 に物質の状態を把握しようとする計算方法である．その結果より，物質の状態量や熱力学的な物性値，

拡散係数などの輸送物性値を評価，検討することができる．実際の観察が困難であるような原子スケ ールでの現象や高速で起こる現象などを把握したい場合には，分子動力学法によるアプローチは適し ており，現在，さまざまな分野において利用されている計算手法である．

分子動力学法もいくつかの種類に分類できる．一つは，古典分子動力学法¹である．この方法では，

計算対象の原子を質点として取り扱い，その原子の運動方程式には，以下の式(3.2.1-1)で示されるニ ュートンの運動方程式(第2法則)を使用する．この方法では，原子の電子構造を取り扱わないので，

例えば，電子のやり取りが行われるような化学反応の計算は取り扱うことができない．

ma

F = (3.2.1-1)

もう一つの代表的な手法として，第一原理分子動力学法 ²が挙げられる．この方法では，量子力学 に基づいた電子状態を計算し，各原子に作用する原子間力を求める．この力を用いて，運動方程式を 解くことにより，分子の運動が分かる．従って，電子状態の変化と原子の運動を結びつけて，理解す ることが可能である．ただし，この手法は，計算の負荷が大きく，そのために取り扱える原子数が限 定されると言われている．近年では，量子力学的に行われる電子状態の計算の一部に経験的なパラメ ータを導入することによって，その計算を高速化したタイトバインディング分子動力学法という手法 も開発されている．この方法によって，電子状態を考慮しながらも，計算対象の分子数を多くできる．

この方法を使用し，表面での化学反応を取り扱った計算も行われている³．

本研究では，計算する系の規模，計算負荷の面から，古典的分子動力学法を採用した．

この方法による計算は，次のようなステップで行われる．

1．分子の配置の設定

2．分子間に作用する力を計算し，それぞれの分子に働く力を求める

3．それぞれの分子において運動方程式を解く，ニュートンの第2法則，式(3.2.1-1)

4．計算された次の時間における位置に分子を移動させる 2-4を適当な回数，繰り返す

計算において，その分子運動を正確に求めるためには，分子モデル，時間間隔の設定，分子間ポテ ンシャルの設定などが重要となってくる．

3.2.2 計算モデル

本研究における計算モデルの概略を，Fig. 3.2-1に示した．上下に鉄原子で構成された壁面を設け，

その間に液体(本研究では，水)を挟む．その後，上部壁面の位置を強制的に移動させていくことによ って，系を減圧させる．従って，初期状態を大気圧にしておけば，負圧の状態を実現できる．水の内 部には，ヘリウム，または，アルゴンの気体分子を配置させ，これを溶解気体とした．

今回の計算では，液体として水を選択した．以下にその理由を述べる．

一つは，水中で起こるキャビテーションにも溶解気体が関係しているからである．Fig.3.2-2 には，

1，5，15 sの各時間におけるキャビティーの画像を示している．滑り速度は，535 mm/s，荷重は，5 N であった．温度は，295Kであり，この温度での粘度は，約0.001 Pa・sである．また，Fig.3.2-3には，

水中におけるキャビティー長さの時間変化を示したものである．これらの図より，キャビティー長さ は，時間とともに成長し，その成長の傾向が第2章2.2に示された潤滑油における実験結果とよく似 ていることが分かる．また，停止後にも，キャビティーの一部が気泡として残り，この点も潤滑油に おける傾向と一致する．従って，水中で起こるキャビティーの成長も，時間とともに起こる水中の溶 解気体の析出と関係していると考えられる．

二つ目の理由に，水中における溶解気体分子の拡散係数などの物性値が把握されていることを挙げ る．計算によるアプローチでは，計算結果を実際の現象と比較し，その妥当性を調査することが必要 である．本計算では，水中での気体分子の拡散係数を比較，評価することによって，計算結果と実際 の現象の妥当性を得た．

三つ目の理由は，水分子にすることで，実験に使用した潤滑油 PAO のような分子量の大きな分子 構造を取り扱う必要がなく，計算負荷が小さくなることを挙げる．この点での有利性も重要である．

3.2.3 ポテンシャル

(1) 水分子

水分子のモデルは，いくつか存在する．ST2，TIP4P，SPCE，CCモデルなどが挙げられる⁴． 本計算では，SPCE モデル ⁵を使用した．このモデルは，Fig.3.2-4 に示されるような構造である．

OH間距離は1Å，2つのOHの結合線の角度∠HOHは104.5°である．O原子の位置に，-0.8476e，H 原子の位置に，+0.4238eの点電荷を持たせている(eは電気素量を示す)．また，この構造において，O，

H 原子の各位置は固定していると考える，すなわち，角度∠HOH は一定である．従って，分子内で のねじれや変角，伸縮運動はしないと仮定している(このような分子モデルを剛体分子モデルという)．

分子間ポテンシャルは，次の式(3.2.3.1-1)のように示される．

∑ ∑ ⁺ ∑∑ ⁺ ∑∑

 +

 





 



 

 



− 

 

 



= 

−

OO O O HO

O H HH

H H OO

O OO

O O

H O

H

r

q q r

q q r

q q r

r

6 12

2

4

2

σ ε σ

φ

(3.2.3.1-1)

右辺の第一項目はO原子間のLennard-Jonesポテンシャル，第二項目はH-H原子間のクーロンポテ ンシャル，第三項目はO-H原子間のクーロンポテンシャル，第四項目はO-O原子間のクーロンポテ ンシャルであり，その和が水分子間に作用するポテンシャルである．各パラメータは，Table 3.2-1に 示す．

(2) ヘリウム，アルゴン

ヘリウム，および，アルゴンは，以下の式で示されるLennard-Jonesポテンシャル⁶を用いた．各パ ラメータをTable 3.2-1に示す．また，2つの分子のポテンシャルの違いをFig.3.2-5に示す．



 





 



 

 



− 

 

 



= 

−

−

−

6 12

4

He He

He He

He He He

He

r r

σ ε σ

φ

(3.2.3.2-1)



 





 



 

 



− 

 

 



= 

−

−

−

6 12

4

Ar Ar

Ar Ar

Ar Ar Ar

Ar

r r

σ ε σ

φ

(3.2.3.2-2)

(3) 鉄原子

上下壁面に使用する鉄は，α-Feとし，そのポテンシャルには，Johnsonポテンシャル⁷を使用した．

Johnsonポテンシャルは，以下の式で示される．原子間距離rijによって，使用する式を選択する．

) 44 . 3 ( 0

) 44 . 3 0 . 3 ( 547967 1

466892 .

0 ) 066403 .

3 ( 115035 .

1

) 0 . 3 4 . 2 ( 581570 .

1 477871 .

0 ) 115829 .

3 ( 639230 .

0

) 4 . 2 9 . 1 ( 436448 .

7 704060 .

2 ) 097910 .

3 ( 195976 .

2

3 3 3

−

=

=

−

=

− +

−

−

=

−

=

− +

−

−

=

−

=

− +

−

−

=

Å

Å Å Å

ij

ij ij

ij

ij ij

ij

ij ij

ij

r

r .

r r

r r

r

r r

r

φ φ φ φ

(3.2.3.3-1)

(4) 水分子-溶解気体分子，水分子-鉄原子，溶解気体分子-鉄原子

水分子と溶解気体分子であるヘリウム，アルゴン間のポテンシャル，水分子と鉄原子間のポテンシ ャル，溶解気体分子と鉄原子⁸間ポテンシャルは，Lorentz-Berthelot の方法 ¹を使用した．この方法，

式(3.2.3.4-1，3.2.3.4-2)により，Lennard-Jonesの各パラメータを決めた．

(

_{i j}

)

ij

σ σ

σ

= +

2

1 (3.2.3.4-1)

j i

ij

ε ε

ε

= (3.2.3.4-2)

(5) 力の計算

ポテンシャルは，分子，原子間の距離を用いて，これまでに示した式で計算できる．運動方程式を 解く際には，ポテンシャルを各x，y，z方向の力に変換する必要がある．以下の式のように，示され，

各座標で微分することによって得られる．この式では，分子iのx方向に働く力を示している．y，z 方向も同様に示される．

∑

=

−

= ^N

j i

ij x

i dx

F d

1 ,

φ

(3.2.3.5-1)

3.2.4 計算方法 (1) 計算時間間隔

計算の1ステップの時間間隔が大きすぎると，細かい分子の動きが再現されないことになる．また，

1ステップでの移動量が大きくなり，分子が他分子を貫通したり，分子同士が接近しすぎるなど，現 実ではありえない結果を引き起こし，これは，計算においても，発散の原因となる．従って，その分 子の状態にあった分子速度を再現できるような時間間隔の設定が必要である．これによって，正確な 分子間力の計算も可能となり，シミュレートされる分子の運動も正確になると考えられる．

これまでの研究報告を調べると，水の液体における解析では，0.5-2 fsが使用されている^9-11．本計 算では，1 fsを計算時間間隔として採用した．

(2) 温度制御

温度Tは，各分子の運動エネルギーの和を基に，以下のように示される¹．

B N

i i

iv Nk

m

T 3

2 2

1

1



 



=

∑

=

(3.2.4.2-1)

m: 質量，v: 分子の速度，N: 自由度の数(3 Ni-1，Niは原子数)，kB: ボルツマン定数(1.38×10^-23 J/K) 従って，計算中に，分子の速度を総和し，式(3.2.4.2-1)に代入することによって，その時間での系の 温度を確認することができる．

本研究では，温度制御の方法として，いくつか提案されている中から速度スケーリング法^12.13を使 用した．速度スケーリング法は，その時間での温度を，式(3.2.4.2-1)によって求め，その値と設定値が 一致するように，速度の補正を行う方法である．従って，系の運動エネルギーを強制的に設定値に一 致させている．設定値との誤差を示す倍率をβとすると，以下のように示される．

v

v′=

β

(3.2.4.2-2)

∑

=

=

=

_N

i i set B set

v m

T Nk T

T

1 2

β 3

(3.2.4.2-3)

この操作を適当な時間間隔で行うことによって，その時間で温度一定の条件が確保される．

本計算では，設定温度を300 Kとし，任意の時間間隔で速度スケーリング法による補正を行った．

(3) 圧力制御

圧力制御にもいくつかの方法が提案されている^13-15．今回の計算では，以下のような方法により圧 力を制御した．

今回の計算の系では，上下壁面に液体を挟んだ構造であることから，上部を可動壁面，下部を固定 壁面とし，上部壁面の位置を移動させることにより系の容積を変化させることができる．このとき，

温度，液体の質量は変化しないので，気体の状態方程式より，圧力が変化することとなる．ただし，

ここで示している圧力は，系全体の圧力であり，局所的なものではない．

(4) 周期境界条件

分子動力学で取り扱う原子数は，その計算負荷に影響するため，無限数に近い原子数を取り扱うこ とは困難である．一般的に，計算に使用される原子数は，アボガドロ数よりも，非常に少ない数であ る．物理量の計算を行う場合には，少ない原子数で行うと，その結果に問題が起こるため，近似的に 原子数を増加させるように工夫した構造が使用される．その一つとして，周期境界条件が挙げられる．

Fig.3.2-6に示すように，実際に計算する系(セル1，セル寸法L×L)の周囲を取り囲むように，その

系がコピーされている状態を考える．セル1 の計算領域のi で示された分子は，セル 2(コピーセル)

のi’と同じである．セル1において jは左のセルに移動する，これは，セル3においてj’で示され，

計算領域のセル1に入ってくることになる．従って，この系の質量は一定に保たれることになる．

周期境界条件をx，y，z方向の全周囲の面に使用することにより，無限の長さを有する立体の計算 領域を作ることができる．

(5) ポテンシャルカット

(4)に示したような周期境界条件の計算領域で計算を行う場合では，一つの分子に対して無限個の分 子との相互作用を考えなければならなくなる．今回計算に使用する Lennard-Jones ポテンシャルは，

Fig.3.2-5 に示すように，ある距離を超えると，そのポテンシャルは非常に小さな値となる．従って，

この距離以降は，計算を無視することが可能である(このような計算の打ち切りを示す距離をカット

オフ距離という)．一般的には，2.0-3.5σにすることが多い^16,17．

また，クーロン力では，式(3.2.3.1-1)に示されるようにポテンシャルが1/rの関数であるので，距離 に対しての変化が，Lennard-Jonesポテンシャルよりも小さい．一般的には，エワルドの方法¹⁸を用い て，この計算の高速化を行うが，今回は，Wolfらによって示されているカットオフ距離を使用する方 法¹⁹を採用し，カットオフ距離25Åを使用した²⁰．

3.2.5 運動方程式の解法

(1) 鉄原子，溶解気体分子の運動方程式の解法¹⁶

鉄原子，溶解分子の運動方程式は，位置の微分値(dr/dt，d²r/dt²)を微小時間間隔における代数式で近 似する差分近似を用いて解く．従って，微小時間間隔で計算し，全分子を移動させる，というステッ プを繰り返すことにより，分子運動を求める．本研究では，速度ベルレ法を使用した．詳細を，以下 に示す．

ある分子における運動方程式は，以下のニュートンの第二法則で示される．

N dt i

md ⁱ

i 1,2･･･

2 2

=

= r

F (3.2.5.1-1)

時間t+∆tにおける位置r(t+∆t)をテイラー展開すると，以下のようになる．

⋅

⋅

⋅

∆ + +

∆ +

=

∆

+ ^{2 2} ₂( )

2 ) ) (

( )

( dt

t d t dt

t td t t

t r r

r

r (3.2.5.1-2)

(3.2.5.1-1)を用いて，(3.2.5.1-2)を書き直すと，以下のようになる．

) ) ) ((

( ) 2

( ) ( )

( ³

2

t m O

t t t

t t t

t ∆ + ∆

+

∆ +

=

∆

+ F

v r

r (3.2.5.1-3)

次に，上で示した時間t+∆tの場合を参考にし，時間t-∆tにおける位置r(t-∆t)をテイラー展開すると，

以下のように示される．

) ) ) ((

( ) 2

( ) ( )

( ³

2

t m O

t t t

t t t

t−∆ = −∆ + ∆ F + ∆

v r

r (3.2.5.1-4)

(3.2.5.1-3)，(3.2.5.1-4)の和をとると，以下のように示される．

) ) ) ((

) ( ( 2 ) ( )

( ² O t ⁴

m t t t t

t t

t+∆ + −∆ = +∆ F + ∆

r r

r (3.2.5.1-5)

従って，時間t+∆tにおける位置は，以下のように示される．ここで，(∆t)⁴以上の項は，無視する．

m t t t t t t

t ( )

) ( ) ( 2 )

( ² F

r r

r +∆ = − −∆ +∆ (3.2.5.1-6)

また，(3.2.5.1-4)，(3.2.5.1-5)の式の差をとると，以下のように示される．

) ) ((

) ( 2 ) ( )

(t+∆t −r t−∆t = ∆tv t +O ∆t ⁴

r (3.2.5.1-7)

従って，上式より，速度vは，中央差分で近似した以下の式で示される．ここでも，(∆t)⁴以上の項 は，無視する．

{

^{( ) ( )}

}

2 ) 1

( t t t t

t t +∆ − −∆

= ∆ r r

v (3.2.5.1-8)

初期条件の速度と位置を与えて，式(3.2.5.1-6)より，次の時間∆tでの位置を計算し，時間tでの速度 は，式(3.2.5.1-8)より求める．この方法をVerlet法という．しかし，この方法では，得られる位置と速 度の時間が∆tだけずれており，∆tでの位置の差より，速度を計算すると，桁落ちが生じる可能性があ る．従って，同時間での位置と速度を用いるように修正されたVelocity Verlet法の方が，計算精度が 高くなる．以下に，その方法について説明する．

式(3.2.5.1-3)を一階微分すると，以下の式で示される．

dt t d t dt

t td t t

t ( )

2 ) ) (

( ) (

2 F

v v

v +∆ = +∆ +∆ (3.2.5.1-9) ここで，速度vの2階微分の項は，無視する．

式(3.2.5.1-1)より，式(3.2.5.1-9)は，以下のように書き変えられる．

dt d t t

m t t t

t F

F v

v( ) ( ) ( ) 2

∆ 2

∆ + +

=

∆

+ (3.2.5.1-10)

前進差分近似より，dF/dtは，以下のように示される．

{

^{( ) ( )}

}

^{( )}

1 t t t O t

t dt

d +∆ − + ∆

= ∆ F F

F (3.2.5.1-11)

従って，時間t+∆tにおける速度は，式(3.2.5.1-10)，(3.2.5.1-11)より，以下のように示される．

{

( ) ( )

}

) 2 ( )

( t t t

m t t t

t ∆ + +∆

+

=

∆

+ v F F

v (3.2.5.1-12)

式(3.2.5.1-3)，(3.2.5.1-12)を用いて，次の時間における速度と位置を計算していく．ここで，(∆t)³以 上の項は，無視する．具体的には，以下のような計算過程をたどる．

・初期条件として位置r0，速度v0を設定する

・力F0を計算する

・時間t+∆tにおける位置r(t+∆t)を式(3.2.5.1-3)から計算する

・時間t+∆tにおける力F(t+∆t)を計算する

・時間t+∆tにおける速度v(t+∆t)を式(3.2.5.1-12)から，計算する

・時間t+2∆tにおける位置r(t+2∆t)を式(3.2.5.1-3)から計算する

・以下，力の計算，速度の計算と繰り返す．

(2) 水分子の運動方程式の解法^1,4,16,21

水分子は，多原子分子であるので，その運動では，重心位置の並進移動と回転移動の2つを考える 必要がある．重心位置の並進運動は，上に示すVelocity Verlet法を使用し，計算を行った．回転運動 の計算は，以下のようにして行った．

回転運動では，次の重心周りの各運動量LとトルクTで示される運動方程式を解く必要がある．

L T dt =

d (3.2.5.2-1)

I

L= ωωωω (3.2.5.2-2) I: 慣性モーメント，ωωωω: 角運動量

分子の回転運動を表す際に，分子の重心位置を原点とした新しい座標系を導入する．ここで，新し い座標系における各座標軸が慣性主軸になるようにとると，慣性テンソルや主慣性モーメントは時間 に依存せず，一つの力学定数として取り扱うことができる．その結果，運動方程式が最も簡単に記述 される．この座標系の概念は，一般に分子座標系といわれる．水分子を分子座標系に固定すると，

Fig.3.2-7のようになる．重心位置を原点にとり，角HOHの二等分線をz’軸に一致させ，2つのH原

子がy’z’平面にくるようにとってある．また，分子座標系と空間座標系の関係は，Fig.3.2-8のように

示される．

分子座標系で運動方程式(3.2.5.2-1)を表すと，各成分において，以下のようになる．

y x y z z z

x z x y y y

z y z x x x

L I L T I

L

L I L T I

L

L I L T I

L











 −

−

=











 −

−

=



 



 −

−

=

1 1

1 1

1 1

&

&

&

(3.2.5.2-3)

式(3.2.5.2-2)を代入し，上式を書き直すと，以下のようになる．この式は，オイラーの方程式といわ れる．

z y x y x z z

y x z x z y y

x z y z y x x

T I

I I

T I

I I

T I

I I

=

−

−

=

−

−

=

−

−

ω ω ω

ω ω ω

ω ω ω

) (

) (

) (

&

&

&

(3.2.5.2-4)

式(3.2.5.2-4)より，各角速度の一階微分は，以下のように示される．

z

y x y x z z

y

x z x z y y

x

z y z y x x

I I I T

I I I T

I I I T

ω ω ω

ω ω ω

ω ω ω

) (

) (

) (

−

= +

−

= +

−

= +

&

&

&

(3.2.5.2-5)

分子座標系と実験室座標系における座標変換には，いくつかの方法があるが，ここでは，オイラー

(Euler)角を使用する方法と，4元数を使用する方法を説明する．

オイラー角をπ，θ，ψとして示す．それぞれの定義は，Fig.3.2-9に示す．この図で示されるように，

3つの回転の操作を行い，2つの座標系を関連付けている．

空間座標系から分子座標系への変換行列Aは，以下のようになる．

 







 







− +

−

−

−

+

−

=

θ φ

θ φ

θ

θ ψ ψ

φ θ φ

ψ ψ

φ θ φ

ψ

θ ψ ψ

φ θ φ

ψ ψ

φ θ φ

ψ

cos cos

sin sin

sin

sin cos cos

cos cos sin

sin cos

sin cos cos

sin

sin sin sin

cos cos sin

cos sin

sin cos cos

cos A

(3.2.5.2-6)

従って，実験室座標系におけるベクトルrは，分子座標系では，以下のようなr’で示される．

Ar

r′= (3.2.5.2-7)

従って，式(3.2.5.2-7)は回転運動における位置変化を示しているので，実際の座標は，重心位置座標 Rとの和になる．従って，以下のような式で示される．

r A R

r = +^t ′ (3.2.5.2-8)

次に，式(3.2.5.2-3)に示したオイラーの運動方程式を，オイラー角と関連付ける必要がある．そこで，

角速度ωωωωをオイラー角によって表すと，以下のようになる．

















+

− +

=















+















+















=















=

ψ θ φ

ψ θ ψ θ φ

ψ θ ψ θ φ ψ θ φ

ω ω ω

&

&

&

&

&

&

&

&

& cos

sin cos

sin

cos sin

sin 0

0

0 0 0

0

2 2

1A A

A

z y x

ω (3.2.5.2-9)

上式を整理して，オイラー角の一階微分を左辺にすると，以下の式(3.2.5.2-10)のようになる．