1 Viscous Fluid Rope の座屈と折りたたみのダイナミクス宮城工業高等専門学校機械工学科永弘進一郎東北大学理学研究科早川美徳ハチミツのような粘性の大きな流体を、平板の上に糸状に垂らすとき、平板の近くで座屈が生じて、流体がコイル状に折りたたまれる現象（Fluid rope coiling）は日常的にもよく経験される。この現象については、定常状態の振動数を説明する Ribe らの数学モデルが知られているが、我々は、既存のモデルに比べて極めて単純で、かつ、transient な状態も含めた全体のダイナミクスを記述できるモデルを提案した。また、臨界レイノルズ数近傍での、Coilling 状態と非 Coiling 状態の転移において、Fluid Rope の振幅はヒステリシスをともなって、不連続に変化することをみいだした。最近、一定速度で並進する平板の上に流体をを垂らした場合に、条件によって様々な振動パターンが生じることが実験的に見いだされた。我々のモデルは、境界条件を変更するだけで、これらの結果をよく再現する。 1 はじめに層流は、流れのレイノルズ数がある臨界値を超えたとき複雑化し、乱流状態へ転移する。一方、低レイノルズ数領域では、強い粘性が、短い波長を持つ固有モードを抑えることから、複雑な現象は生じにくい。しかし、流体の表面が拘束されておらず、自由に変形できる場合は、例外的に、低レイノルズ数側で生じる不安定性が存在する。ハチミツをトーストの上に細く垂らすと、パンの上で蜜の “ロープ”（以下、ﬂuid rope と呼ぶ）が座屈し、コイル状に折りたたまれる。この現象は “Fluid rope coiling” と名付けられていて、流体がある程度以上大きな粘性を持つときに初めて現れる、低レイノルズ数領域での不安定化現象である。日常的にもよく経験されるこの座屈不安定化現象について、半世紀ほど前から様々な実験や理論的考察が行われてきた。最も古くは、Barnes と Woodcock によるエンジンオイルを用いた実験がある [1, 2]。彼らは、折りたたみの振動数 Ω が、落下高さ（ノズルの出口と平板間の距離）H に比例して増加することを報告した。その後、Cruickshank らは、ﬂuid rope が自発的な回転を起こすための最小の高さ H ∗ が存在し、その付近で振動数が H の減少関数となりうることを確認している [3, 4]。また、Ribe と Mahadevan らはそれぞれ、Fluid rope coiling の定常状態における折りたたみ運動の振動数を、Rope 内に生じる応力の単純なバランスを考えることによって説明した [6, 7]。定常状態では、座屈点近くにおいて、駆動力と流体の粘性応力 σV の釣り合いが実現している。駆動力として、ノズルから流体が押し出される応力 σP 、重力による応力 σG 、慣性による応力 σI がある。流体の落下高さが低い場合、重力の効果は小さいので、σV ∼ σP という釣り合いが実現する。落下高さが徐々に大きくなり、流体の吐き出し口と座屈点が離れると、σP は重要でなくなり、かわりに σV ∼ σG が実現する。吐き出し口がさら高くなると、流体は重力によって引き延ばされ、断面積は小さくなる。重力の効果は、断面積との積で現れるので、σG の効果は小さくなり、替わ 2 りに σV ∼ σI というバランスが実現する。上記それぞれの釣り合いから、折りたたみの振動数が ΩV ∝ QH −1 a−2 , ΩG ∝ g 1/4 (1) 3/4 −2 −1/4 Q a ν , (2) ΩI ∝ Q4/3 a−10/3 ν −1/3 (3) となることを示すことが出来る。ここで、Q : 体積流量、ν : 動粘性係数、g : 重力加速度、H : 落下高さ、a : 座屈点での Rope の半径である。これら 3 つのスケーリング領域が確かに存在することは、実験的にも確かめられている [8, 9]。これらの事実にたいして、Fluid rope の理論的な扱いは、線形安定性から H ∗ を説明する方法 [4, 5] や、定常状態に限定して振動数を記述する数値理論が知られている [7]。しかし、垂直に流れ落ちる Fluid rope が座屈し、定常的な折りたたみ運動へと移るプロセス再現し、座屈の臨界 Reynolds 数まで説明する理論はまだ無い。さらに最近の実験から、一定速度で並進する平板の上に Fluid rope を落とすと、様々な振動パターンが現れることが報告された。並進速度と落下高さをパラメータに選ぶと、10 種類以上に分類される振動状態が観測される。本稿では、上記のような複雑で多様な現象の全体を記述可能な、単純な微分方程式モデルを提案し、実験事実の説明を試みる。 2 モデル方程式の導出以下では、ﬂuid rope の自発的な回転を記述する擬一次元的なモデル方程式を導出を行う。図 1 に、ノズルから落下する流体の様子を示す。ノズルから d 落下する流体は、下部の平板に衝突して急速に減速する。モデルでは、この減速と圧力勾配によって生じる軸方向の粘性応力 σt が、流体を座屈させ回転運動を起こすもっとも重要な駆動力であると見なし、表面張力や粘弾性の効果を考えない。その上 z S で、数学的な記述を単純化するために、以下の２点を仮定する。 x 1. ﬂuid rope の変形は小さいものとし、２次の微小量ヒネ y リと曲げについてのスティフネス*1 を無視する。 2. 折りたたまれた ﬂuid rope が折り重なっている部分（図１参照）の最上部は、供給される流量によって上下動する。しかしここでは、最上部の上向きの成長速度は、そ図1 の部分の流体の下向きの流速と常に釣り合っていて、静性流体の折りたたみの様子ノズルから落下する粘止しているものとする。以下、仮定 2. を元にして、座標の原点を図 1 に示したようにコイル部分の最上部におき、この原点は動かないものとする。z 軸をノズルと平行にすえ、反対方向を正にとる。以後、z > 0 での流体の運動のみを考える。仮定 1. より、ﬂuid rope の接線ベクトル *1 ここでは、曲げスティフネスとは、ﬂuid rope をある有限の速度で曲げるときに、中心線の外側が引き延ばされ、内側が圧縮されることによって生じる粘性応力を指す。 3 t は、z 軸にほぼ並行であるので、t 方向の微分を ∂ に置き換える。rope の断面積を S 、流速を ∂z w として、連続の式は ∂S ∂ = (Sw), ∂t ∂z (4) と書ける。流体の密度を ρ として w の運動方程式は ( ∂ ∂ +w ∂t ∂z ) w= 1 ∂ (Sσt ) − g, ρS ∂z (5) である。この二つの方程式は、ﬂuid rope が完全に直線である時の流れを記述している。粘性によって座屈が起こると、rope の中心軸は z 軸からずれる。その変位を q = (qx , qy )、速度を u = (ux , uy ) とすれば、両者は ( ∂ ∂ +w ∂t ∂z ) q i = ui , (6) という関係で結ばれる（なお、i = x, y ）。図 1 に抜き出して示してある ﬂuid rope の微小要素に加わる力は、σt の xy 平面への射影成分と、隣り合う微小要素間にはたらく剪断応力 σi である。粘性 ∂(ui + wti ) と書けるから、ui の運動方程式は次の形に書ける。 ∂z ( ) ∂ ∂ 1 ∂ +w ui = {S (σi + ti σt )} . (7) ∂t ∂z ρS ∂z 係数を η として、剪断応力は σi = 軸方向応力 σt は、 σt = −p − 2η ∂w ∂z という形を持つが、圧力 p は、rope 表面での半径方向の応力が σrr = 0 である条件から、以下に ∂ur である。体 ∂r ∂w ∂w ∂w ∂ur =− から、p = η である。よって、σt = −3η となる [15]。積保存の条件 2 ∂r ∂z ∂z ∂z ノズルの位置 z = H での境界条件は、 π qi (H) = 0, ui (H) = 0, S(H) = , w(H) = −|win | 4 示すように消去できる: ﬂuid rope の半径方向の流速を ur として、σrr = −p − 2η とし、z = 0 にて qi′ (0) = 0, u′i (0) = 0 (8) を課す（プライムは z 方向の微分を表す）。流速 w(z) の z = 0 での境界条件は、有限の値を与えなければならないが、これを指定するために、現象論的なパラメータを導入する。ノズルから落下する流体は、はじめ引き延ばされ、平板付近では圧縮される。Cruickshank らは、伸張から圧縮へとうつる高さ z = ξ において座屈が生じると考え、流速 w(ξ) を実験的に測定した。その結果、 w(0) = βw(ξ), (β = 0.14) (9) という関係が、その他の流れの条件に殆どよらずに成り立つことをみいだしている [16, 17]。この関係を w に対する境界条件として用いる。方程式 (4)、(5)、(6) および (7) を無次元化して得られるレイノルズ数 Re = d|win |ν 、フルード 2 数 Fr = win /gd、落下高さ H/d をモデルのコントロールパラメターとする。本稿にて以下に示す結果は、空間刻みを ∆z = 0.1、時間刻みを ∆t = 0.25Re∆z 3 としてオイラー差分を用いて計算した。初期条件は、w(z) = −1、qi (z) = ui (z) = 0 とするが、qi (z) には、振幅 0.01 程度のラフネスを与えている。 4 (a) (c) coiling radius 4 10 3 (b) 10 Fr 0.01 qy (0) 1 2 10 6 3 0 1 2 Re 3 0 0 1 10 Coil -1 -0.01 -1 0 1 -0.01 qx (0) 0.00 0.01 1 Uncoil H/d=10 1 2 qx (0) Re 5 10 図2 (a) 座屈が生じ、回転運動が安定な条件での Rope 最底部の軌道 (Re = 1, Fr = 1, H/d = 10)。(b) 座屈が生じない場合の軌道 (Re = 3, F r = 1, H/d = 10)。(c) 座屈が起こる領域と起こらない領域の状態図。プロットの△と▼はそれぞれ、レイノルズ数の上昇時と下降時の臨界値を表す。 3 座屈の生じる条件 Fluid rope の座屈は、レイノルズ数がある臨界値 Re∗ を下回ったときに起こるが、この臨界値は F r と H/d に殆ど依存せず Re∗ ≅ 1.2 程度であることが実験的に見いだされている [3]。モデルの数値解から得られた Fluid rope 底部の軌道 (qx (0), qy (0)) を図 2(a)(b) に示す。Re = 1 の場合、軌道は振幅が１程度の円形の軌道に収束するが、Re = 3 では、軌道は原点へと収束している（座屈はおこらない）。実験との比較のため様々な F r の値について、臨界レイノルズ数をプロットしたものが図 2(c) である。1 < Fr < 104 の範囲にわたって、Re∗ は Fr に依存して変化するが、変動は一桁に収まる程度に弱く、値も実験とほぼ一致している。また、Fr > 100 では、臨界レイノルズ数に履歴効果が現れる。図 2(c) 内に挿入したグラフが示すように、Re の関数としての Fluid rope の振幅は、Re が上昇するときと下降するときで、不連続店が異なる。つまり、Coiling 状態への転移が、subcritical 分岐として起こることを示唆している。 4 定常状態での折りたたみ振動数この節では、臨界 Reynolds 数の定常状態における折りたたみ運動の振動数について議論する。重力の効果が大きい場合 (Fr ≪ 1) と小さい場合 (Fr ≫ 1) に分けて考える。 4.1 モデルによる結果図 3 に、定常状態におけるコイリングの振動数 Ω と落下高さ H の関係をしめす。(a) 重力が小さい場合 (Fr > 100)、H ≅ 5 を境に、振動数は H −1 に比例する減少から、Ω ∝ H に従う増加へ転じる。(b) 重力が大きい場合 (Fr < 0.05) は、振動数は、H > 3 の範囲で単調増加的となり、振動数は H 2 に比例して増加する。これらの結果は、Ribe らによる実験結果（挿入図）とかなりの部分で一致している。しかし重力が強い場合の Ω の指数は実験と比べて小さい値であり、また 5 10 8 6 10 3 5 (s-1) 10 2. (s-1) 2 10 2 d = 0.5 cm 0.5 1 2 5 1 10 1 H 20 2. 4 0.2 0 d/g = 0.02s 50 -1 (a) Re=3 10 1 (b) Re=0.1 3 100 10 30 H/d H/d 図 3 コイリング振動数 Ω と落下高さ H/d の関係（挿入図は Ribe らによる実験の結果 [9]）。(a) Re = 3, Fr = 100(⃝), 200(△), 400( )。(b) Re = 0.1, Fr = 0.01(⃝), 0.02(△), 0.05( ) H = 7 程度で生じている振動数のジャンプも再現されない。 4.2 コイリングの振動数とスケーリング則以上の結果を、三種のスケーリング則 (1)、(2)、(3) と比較するためには、座屈点での Rope の半径 a を決定しなければならない。ここでは、ロープが伸張から圧縮へと移る高さ z = ξ と座屈点が等しいとする Cruickshank の考えを採用して a を見積もる。z = ξ では流速 w(z) が極大になっているはずであるから、この点を z 軸の原点に取り直し、境界条件を w′ (0) = 0, w(H) = −1 とする。その上で無次元化した (5) 式における定常状態 dw 3 d w = w dz Re dz ( 1 dw w dz ) − 1 . Fr (10) を考える。1/Fr = 0 ならば、この式は w(z) = −1 という自明解をもつので、Fr ≫ 1 の場合を考え、摂動解 w(z) = −1 + 1 φ(z). Fr を求める。これを (10) に代入して得られる微分方程式は φ(z) について簡単に解くことが出来、 w(z) = −1 + } 1 { z − H + (e−z − e−H ) . Fr を得る。定常の連続の式 Q = −Sw(z) から、半径 a は √ a= 1 w(H) =√ . w(0) 1 − (1 − H − e−H ) /Fr (11) となる。まず、Fr ≫ 1 である場合から考える。もし H ∼ 1 程度の時は、Ω ∝ ΩV が成立する。(11) 式から a ≅ 1 であるので、Ω ∝ H −1 を得る。落下高さが H ∼ Fr 程度に高くなると、流れに重力の効果が無視できなくなる。このとき Ω ∝ ΩG が成り立つとすれば、(11) 式は a−2 ∼ 1 + H/Fr を与える。よって、Ω ∝ H + const. という、Barnes らの実験によって知られている関係が再現される。 Fr ≪ 1 では、座屈の生じるレイノルズ数もまた小さい。よって (10) 左辺の慣性項は無視でき、鉛直方向の流れは、重力と粘性のバランスから決定される。H が高くなるにつれ、Rope は強く引 6 き延ばされ、その振る舞いは漸近的に a ≅ (Qν/g)1/2 H −1 であると導くことが出来る [7]。この関係を Ω ∝ ΩG と Ω ∝ ΩI にそれぞれ代入してみると ΩG ∝ H 2 , ΩI ∝ H 10/3 (12) と類似した指数を持つ関係が得られる。我々のモデルの数値解から得られた指数は 2.0、Ribe らの実験では 2.5 が観測された。Ribe らは、これらの指数が彼らの数値計算とも一致することから、スケーリング領域 Ω ∝ ΩI が確かに存在することを主張している。これらの指数が (12) 式のどちらを反映したものなのかを厳密に判定するには、振動数の重力依存性を見ればよい。ΩG は重力に依存し、ΩI は依存しない。a ≅ (Qν/g)1/2 H −1 より ΩG ∝ g 5/4 である。我々のモデルの結果は、 10 < Fr−1 < 103 の範囲において Ω ∝ Fr−1.25 の関係をあたえる。よって Fr ≪ 1 の条件の下、我々のモデルはスケーリング則 Ω ∝ ΩG を再現していることが分かる。落下する流体の慣性と粘性応力の釣り合いによって説明されるスケーリング領域 Ω ∝ ΩI や、実験で確認されている振動数の不連続な飛びを、我々のモデルは再現しない。実験において Ω ∝ ΩI が観測される領域では、ノズルから落ちる流体はほそく引き延ばされていて、 102 Hz から 103 Hz 程度の非常に高い周波数で折りたたまれている。このような運動状態では、我々が (9) 式で仮定したような z = 0 での境界条件は成り立っておらず、Rope の座屈点は激しく移動し、β の値も一定にはならないだろう。このことが、我々のモデルが高い振動数の領域を再現しない要因の一つとなっていると考えるが、座屈点の振動をモデルの境界条件へ効果的に取り込むことは、難しい課題である。 5 並進ベルト上での Meandering 不安定性 Webstar と Lister は、速度 U0 で水平方向に並進する平板の上での Fluid rope の座屈を実験的に観測した [11]。定常なコイリング状態にある Fluid Rope について、平板を小さな速度で並進させる場合は、 Rope の運動は大きな影響を受けず、回転を続ける。よって、平板上にはサイクロイド状の「わだち」がのこる。これを Translated coiling(TC) 状態と呼ぶ。平板を徐々に加速すると、平板の運動方向に平行な運動が突如として消え、「わだち」はサインカーブ状の曲線となる [Meandering(M) 状態]。U0 が十分に大きいと、もはや振動は消えてしまい、Rope は静的な懸垂線を形成する。また、落下高さが H/d > 7 の領域では、Rope は複雑な振動を呈し、「8 の字型」や「W 字型」の曲線が現れる。ある条件においては、平板の並進方向に対して左右対称性の破れた振動が現れることも分かっており、興味深い。我々のモデルにおいて、平板が並進している状態を再現するために、ux (z) の下端での境界条件 (8) を次のように変更する。 ux (0) = γ (13) ここで γ は、Rope が平板の並進方向へ引かれることに起因する速度勾配であり、γ > 0 とする。実験において調節可能な量は、平板の並進速度である。境界条件 (13) は、平板上のある高さにおける並進方向の速度勾配を固定することを意味する。よって、実験との比較を定量的に行うためには、γ と平板の速度の関係が明らかでなければならないが、残念ながらそれを求めることは出来ていない。ここでは、γ が U0 の一次に比例すると考えて、議論を進める。以後、文献 [11] での実験 7 H=10 0.8 (b) (a) (c) 0.4 y 0.0 -0.4 -0.8 -1.0 s=5 s=10 -0.5 x 0.0 -1.0 s=6 -0.5 0.0 -1.0 -0.5 0.0 x x 図 4 定常状態における Rope 下端の軌道 (H = 10)。(a) Translated coiling (γ = 4). (b) Meandering (γ = 10). (c) H > 7 の領域で現れる Translated coiling と Meandering の混合状態。 40 10 (a) 9 (b) 8 Belt speed (cm/s) 30 Catenary 20 M Crossover 10 7 6 5 catenary meander figure 8 trans. coiling double meander disordered stretched coiling slanted loops W state double coiling 4 3 2 TC 0 5 10 15 H/d 図5 1 20 0 2 3 4 5 6 Nozzle height (cm) 7 8 並進平板上に落ちる Fluid rope の振動状態図。(a) モデルによる結果。(b) 実験的に得られた状態図 (d=0.8cm) 文献 [13] より転載。に対応させ、νQ/gd4 = 3.3 × 10−2 、gd3 /ν 2 = 3.3 × 10−3 とする。 γ = 0 で Rope が定常な回転運動をする場合を考える。γ が小さければ、図 4(a) に示すような、原点を内側に含むわずかにひずんだ円形の軌道が現れ、これは TC 状態と言える。γ がある値を超えると、M 状態に対応する細い 8 字型の軌道が、x の負の方向へ後退した位置に実現する [図 4(b)]。また、落下高さが H/d = 8 を超えると、上記 (a) と (b) のクロスオーバーに対応する軌道が現れ、H が大きくなるに従い、その領域は拡大する。この三つの状態を H-γ 平面にプロットすると、実験的と類似した状態図が得られる。実験では、H ≅ 6.5cm, U0 ≅ 3cm/s におおいて、M 状態と TC 状態の両方が観測される領域が現れ始め、それよりも高い H で、複雑な運動状態が観測される。このことから、並進平板上での Fluid rope の複雑な運動状態は、M 状態と TC 状態の混合によって起こっていると予想できる。しかしこの描像から、左右対称性の破れた振動状態などを具体的に説明することは難しく、今後の課題として残る。 6 まとめ粘性流体の座屈と折りたたみという、半世紀近く研究の対象となっているよく知られた現象に対して、もっとも単純なモデルを導出し、それが実験事実をよく再現すること示した。さらに、最近の実験によって報告された、並進する平板上での Meandering についても、現モデルの境界条件を 8 変更するのみで、その振る舞いが説明できた。参考文献 [1] G. Barnes and R. Woodcock, Am. J. Phys. 26, 205 (1958) [2] G. Barnes and J. MacKenzie, Am. J. Phys. 27, 112 (1959) [3] J. O. Cruickshank and B. R. Munson, J. Fluid Mech. 113, 221 (1981) [4] J. O. Cruickshank, J. Fluid Mech. 193, 111 (1988) [5] B. Tchavdarov, A. L. Yarin and S. Radev, J. Fluid Mech. 253, 593 (1993) [6] L. Mahadevan, W. S. Ryu and A. D. T. Samuel. Nature (London), 392, 140 (1998); 403, 502 (2000) [7] N. M. Ribe, Proc. R. Soc. Lond. A 460, 3223 (2004) [8] N. M. Ribe, M. Habibi and D. Bonn, Phys. Fluids, 18, 084102 (2006) [9] M. Maleki, M. Habibi, R. Golestanian, N. M. Ribe and Daniel Bonn, Phys. Rev. Lett. 93, 214502 (2004) [10] M. Habibi, M. Maleki, R. Golestanian, N. M. Ribe and D. Bonn, Phys. Rev. E 74, 066306 (2006) [11] S. Chiu-webstar and J. R. Lister, J. Fluid Mech. 569, 89 (2006) [12] N. M. Ribe, J. R. Lister and S. Chiu-webstar, Phys. Fluids, 18, 124105 (2006) [13] S. W. Morris, J. H. P Dawes, N. M. Ribe and J. R. Lister, Phys. Rev. E 77, 066218 (2008) [14] S. Nagahiro and Y. Hayakawa, Phys. Rev. E 78, 025302(R) (2008) [15] F. T. Trouton, Proc. R. Soc. Lond. A 77, 426 (1906) [16] J. O. Cruickshank and B. R. Munson, Phys. Fluids 25, 1935 (1982) [17] J. O. Cruickshank and B. R. Munson, Phys. Fluids 26, 928 (1983) [18] J. Teichman and L. Mahadevan, J. Fluid Mech. 478, 71 (2003)