数理解析研究所講究録 1413 巻 2005 年 130-138130

膜の座屈による皺パターン
Wrinkle Patterns on Buckled Membranes
内田就也

Nariya UCHIDA

Department of Physics, Tohoku University

1

緒論

弾性体の座屈問題は Euler の安定性解析に始まる長い歴史を持つ [1]. その最も基本的な例である, 軸に平行な圧縮力を受けた細い棒は, 力が長さの逆 2 乗, 直径の逆4 乗に比例する臨界値を超えたとき, 曲けに対して不安定となる. 座屈したあとの棒の平衡形状は, Euler による解析解が elastica として知られている. 棒に比べると 2次元弾性体, すなわち膜の座屈形状の解析は少し 困難である. 膜の大変形を記述する方程式は 4 階非線形連立偏微分方程式となり, その解は数値的にのみ得られることが多い. 一方, 我々の身のまわりにおいて, 弾性膜の座屈問題は, 乗物のボディや耐圧容器の設計に際して重要となるほ力\searrow 紙をくしゃくしゃに丸めたときの皺のバターンにも関連がある. また皮膚の皺なとも, 基板上に接着した膜の座屈として捉えることができる. 近年, このような皺の力学に対する物理からのアプローチが盛んになって来ている [2, 3, 4, 5, 6]. 力学的に見たとき, 弾性膜の大変形パターンを構成する基本要素は, 弾性エネルギーが集中する線状領域 (リッジ) である. リッジは紙の皺における折れ目 (山折れ, 谷折れ) に相当する. ただし紙の折れ目は塑性変形して鋭角的になっているのに対し , 弾性膜のリッジは曲け弾性のために丸みを帯ひ,その曲率は, 曲けと面内歪みの 2 種類の弾性エネルギーのバランスで決まる. リッジ(尾根) という呼称は, それが再端点をビークとする吊尾根のような形状を持つことによる. Witten ら [2] のスケーリング理論によると, 長さ のリッジの弾性エネル
ギーはに比例する. 多数の同様なリッジを含む皺パターンでは, 1 本のリッジあたりの面積はだから, 単位面積あたりの弾性エネルギーはに比例する.
本論文では, 弾性膜の座屈によって生じる皺パターンの特徴を議論する. これはリッジ間の相互作用の解析に帰着するが, 特にここで着目するのはリッジの方向自由度である. 四囲から等方的な圧力を受けて座屈する膜においては, 皺が特に走りやす$\langle$$\text{ー}$$\ell$$\ell^{1/\}$$\ell^{2}$$\ell-5/3$

131

い方向はな , すべての方向が等価である. 自発的対称性の破れによって特定の方向に皺が生じ, 他の皺と合体成長して座屈が進行する. この過程はベクトルスビン系の相転移ダイナミクスとのアナロジーで理解することができる. ここでスビン変数と$\langle$なるのは膜の勾配ベクトルであり, シアー歪みはスピン間の有効長距離相互作用を媒介することが示される. 本論文ではシアー歪みの役割を明示するため, 弾性膜と ,シアー弾性を持たない模型的な「流体膜」を比較しながら議論を進める.

2

モアル

本論文では弾性膜の可逆変形を, 膜厚が変形の特徴的な長さに比べて十分小さい薄膜極限で取り扱う. 変形する前の膜は x 架平面上にあったものとして, 膜を構成する各点の変位を $U=U(x,y)=(u, h)=(\mathrm{u}_{t}, u_{\mathrm{y}}, h)$ とする. また, 変形によって生1 と仮定すると, 力学的平衡条件下においじる膜の勾配は小さい, すなわち |て歪みの水平成分
$h|$

(1)

$U_{1j}.= \frac{1}{2}(\partial_{}u_{j}+\partial_{j:}u)$は $O((\nabla h)^{2})$ の微小量になることが, 三平方の定理から容易に導かれる. このことを利用して Cauchy の歪みテンソルを水平成分について線形化し,(2)$(i,j=x,y)$$E_{1\mathrm{j}}.= \frac{1}{2}[\partial_{1}.u_{jj:}+\partial u+(\partial_{}h)(\partial_{j}h)]$) 近似である [7]. この近似でとするのが, いわゆる-von K\’am\’an (以下の和として表と曲けエネルギーは膜の弾性エネルギー $F$ は歪みエネルギーされ, 前者は$\mathrm{F}\mathrm{v}\mathrm{K}$$\mathrm{F}\tilde{\mathrm{o}}\mathrm{p}\mathrm{p}\mathrm{l}$$F_{l}$$F_{b}$$F_{*}= \int dxdy(\frac{\lambda}{2}E_{\ddot{u}}^{2}+\mu E_{1j}^{2}.)$(3)の形を取る. ここに と は 2 次元 Lame’ 定数であり, $K=\lambda+\mu$ が面内での等方的な変形 (面積変化) に対する弾性率, が面内シアー弾性率である. 曲け弾性エネ
ルギーは膜の平均曲率 H=(\nabla 2h)[1+O( h)2] の 2 乗に比例するが, 再ひ $O((\nabla h)^{2})$$\lambda$$\mu$$\mu$の項を無視するならば$F_{b}= \frac{\kappa}{2}\int dxdy(\nabla^{2}h)^{2}$(4)で与えられる. 曲け弾性率 は膜厚の 3 乗に比例することが示されている.
次に変形のダイナミクスについて, 本論文では o rdamped な方程式$\kappa$$\frac{\partial \mathrm{u}}{\partial t}=-L_{\mathrm{u}}\frac{\delta F}{\delta u}$,$\frac{\partial h}{\partial t}=-L_{\hslash}\frac{\delta F}{\delta h}$

(5)

132

(a)

(b)
$\mathrm{h}*\mathrm{M}$$\mathrm{t}-\{\infty$$\iota-u\mathrm{o}0$g’ 一一

G 龜 man \mbox{\boldmath $\omega$}\tilde ml 騙 n$1\approx 400$$\mathrm{t}-\mathrm{t}00$$\mathrm{t}\cdot \mathrm{O}\mathrm{O}$t-160O$\mathrm{w}/u\mathrm{e}$$\mathrm{h}\cdot[] \mathfrak{g}[][][]$$1\cdot 6\cdot 00$$\mathrm{t}\sim 1\mathfrak{W}$$.\mathfrak{g}[] u\mathrm{e}\mathfrak{g}r\mathrm{a}\mathrm{r}m$$\mathrm{t}[]\alpha 0$

図 1: 膜の高さ$h$

伽wJan cur 鴨加糟

$\mathrm{b}6\cdot 00$のスナップショット (時間発展) , およひ時刻 $t=6400$ における 2 乗勾配と Gauss 曲率 $\det$ ( き $h$ ). (a) 弾性膜 $(K=2, \mu=1)$.

$(\nabla h)^{\}$.

流体膜 $(K=2, \mu=0)$ グレースケール
は, 高さと Gauss 曲率は正 (負) の所が明色 (暗色), 2 乗勾配は 0 の点が黒となるよう取った.
$(\mathrm{b})$を仮定する. このような模型的なダイナミクスを採用する理由には, 数値計算の簡単化のためという側面もあるが, 後述するように, 相転移ダイナミクスとのアナロジーからスケーリング則を検証することができるという利点もある. 制御パラメータとしては水平歪みの空間平均 . を取る. 時刻 $t=0$ に等方的な圧縮 $U\text{。}=U\delta_{-j}(U<0)$$F_{1j}$を与えて固定し, その後の座屈パターンの時間発展を見ることにする. これは弾性膜を枠にはめて, 瞬間的に枠を狭めることに対応している. シミュレーションにおいては $U(x,y)$ に周期境界条件を仮定した. 用いたパラメータ値は, 格子長, 時間ステップをそれそれ長さ , 時間の単位に取ったとき, 格子サイズ, 時定数 $L_{h}=L_{\mathrm{u}}=0.01$ ,$\kappa=0.01$ である. シアー弾性の役割を調べるため, 弾性膜 $(\lambda,\mu)=(1,1)$ と『流体膜」 $(\lambda,\mu)=(2,0)$ の 2 つの場合についてシミュレーションを行なった.
$512^{2}$

133

(b)

$u’$

$\mathrm{u}’\sim$.

図 2: (a) 曲け弾性エネルギー及ひ (b) 高さの 2 乗平均の時間発展. 実線はべき乗則によるフイツテイング. fluid (fast u) とあるのは流体膜に対する有効的な動力学方程式 (10) を解いた結果.

3 成長ダイナミクスを図 1
弾性膜と流体膜のそれそれについて, 座屈後の皺パターンの時間発展ののプロットから, 弾性膜では線状の(上段と中段) に示す. 図 1 下段に示す
リッジ, 流体膜では等方的な丘状の構造がパターンの基本単位であることが分かる.
また弾性膜では Gauss 曲率のゆらきが抑えられていることに注意する. 皺の成長の$\mathrm{g}.\neq$$(\nabla h)^{2}$に緩和して歪みエネルギーを解放すが平衡値初期においては, 勾配 lる一方, 曲けエネルギーは増加する. その後, 曲けエネルギーを減少させるため$-パ
およひ高さ の 2 乗平均を図 2 に示
ターンは粗大化していぐ曲けエネルギー$h|$$\sim\sqrt{|U|}$$h$$F_{b}$す ( グラフから分かるように, これらの量は粗大化過程を通じてべき乗則(6)$h\propto t^{\beta}$$F_{b}=t^{-a}$てよくフィットされる. べき指数は弾性膜では $(\alpha, \beta)=(0.47,0.54)$ , 流体膜では$(\alpha,\beta)=(0.47,0.45)$ であった.

弾性膜に対するべき乗則は, 前述のリツジ描像によって説明することができる. 実
際, 弾性エネルギー $F$ が平均のリッジ長 の 3/5 乗に比例することと, overdamwd, したがってなダイナミクス $dL/dt\alpha F/L$ とを組み合わせれば,結果はス$6/11=0.545\ldots,$ $\beta=5/11=0.454\ldots$ となる [3]. 上の数値シミュレーシ 1$L$$L\propto t^{3/11}$$\alpha=$$\sqrt[\backslash ]{}$

ケーリング描像を支持していると言える.
次に流体膜のべき乗則を理解するために, 水平変位 は垂直変位眉こ隷属する速い変数であると仮定して, これを断熱的に消去する. 具体的な計算は次節に示すが,その結果は勾配ベクトル$u$$m=\nabla h$(7)134に対する有効的な運動方程式として次の形に表すことができる.
$\frac{\partial m}{\partial t}=\Gamma(\overline{U_{i*}.}+\frac{1}{2}\overline{m}^{\mathrm{F}})\nabla^{2}m-\Gamma\kappa\nabla^{2}\nabla^{2}m$

この式は, 有効$J\backslash$.

(8)

ミルトニアン(9)$H=H[m]= \frac{1}{2}\int dr[\lambda(\overline{U_{i}}+\frac{1}{2}\overline{m}^{T})m^{2}+\kappa(\nabla m)^{2}]$の形に書き直すこともできる. これは相転移ダイを用いて,ナミクスで言うところの保存系モデルの式とよく似ているが, ハミルトニアンの 4 次の非線形項が空間平均を含んでいる点が大きく異なる. 実際, このグローバル結合のため, この系にはベクトルスビン系のような位相欠陥が存在しない。粗大化過程の後期において, 方程式 (9) の右辺第一項は空間各点でほぼ同じように 0 に近付く. このことから次元解析が適用できて特微的な長さが $L\sim t^{1/4}$ となることが予想される. これは $\alpha=\beta=1/2$ を意味するが, 数値的に得られた指数は 1/2 より若干小さい. このすれの理由は明らかではないが, 2 吹元 Laplace 方程式の基本解が距離に対数的に依存することと関連した有限サイズ効果である可能性がある. なお,上の有効モデルを直接数値的に解いた結果は図 2 において fluid (fast u) として表示している. これは, 成長過程の初期を除いて, 最初のモデルの結果とよく一致しており, 水平変位を速い変数として消去することを正当化している.
$\partial m/\theta t=\Gamma\nabla^{2}\delta H/\delta m$$\mathrm{X}\mathrm{Y}$

4

方向相閤

本節では, 水平変位 を断熱的に消去した有効モデルを用いて弾性膜の皺バター
自由エネルギーは について線形なので, 力学的つりあいの方ンを解析する.
程式 $\delta F/\delta u=0$ は について完全に解くことができ, その結果を $F$ の表式に代入することで, 有効自由エネルギーが膜の勾配 $m$ の汎関数として得られる. その結果, 歪みエネルギーは次の形になる.
$u$

$u$

$\mathrm{F}\mathrm{v}\mathrm{K}$$u$$F_{\epsilon}$$= \int dr[\frac{\lambda}{2}\{$ $( \sigma_{||}^{-=}..+\frac{1}{4}m)m^{2}\}+\mu\{$ $( \sigma_{1j}^{-}.+\frac{1}{4}\varpi_{1}.7\Gamma_{j})m:m_{j}\}$$+ \frac{\mu(\lambda+\mu)}{2(\lambda+2\mu)}\{m^{2}-\frac{\nabla_{1}\nabla_{j}}{\nabla^{2}}(m:mj)\}^{2}]+\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{s}\mathrm{t}$

.

(10)

ここで被積分因子の第 1 行目は前節の流体膜の場合と同じ, 非線形項がグローバルは Laplacian結合となったモデルの形をしている. 一方第 2 行目においての逆演算子であるが, 膜の Gauss 曲率 det( き ) を用いてdet( き h) と書き直せることに注意する. すなわち第 2 行目は, 膜の各点でのGauss 曲率の間に, シアー歪みを介した長距離斥力が働くことを意味する. このことから, 弾性膜において Gauss 曲率のゆらきが抑制されるという前節のシミュレーシ a$1/\nabla^{2}$$\mathrm{X}\mathrm{Y}$$h$$(\eta\nabla^{2})$$m^{2}-\nabla_{1}.\nabla_{j}/\nabla^{2}(mtmj)=$135(b)(a)におけるスナ プショッ.
図 3: (a) 弾性膜およひ (b) 流体膜の相対方向相関関数
トから計算したもの (10 サンブルの平均). 縦軸は G。I(0) $=1$ となるよう規格化 てある.
$G_{\mathrm{r}\epsilon 1}(X,Y)$$t=1\alpha \mathrm{n}$$\nu$$\mathfrak{U}$

ン結果が理解される. この有効相互作用はまた, リツジの方向相関に影響を与える.
これを示すためにます, リッジの方向は概ね $m$ に垂直であり, また, 反転 $marrow-m$に対して不変であることに注意する. このことから, 方向自由度をテンソル(11).

$Q_{jj}=m:m_{j}- \frac{1}{2}\delta_{1j}$の形て取り出すのが都合がよい. これは $m=m(\cos\theta,\sin\theta)$ とすれば $aQ_{eae}=-Q_{w}=$$(m^{2}/2)\cos 2\theta,$ $aQ_{eu}=Q_{yx}=(m^{2}/2)\sin 2\theta$ で与えられる, 対称かつトレースが零のテに比例すンソルである. 有効自由エネルギー (11) を Fourier 表示すると,については,る項が現れるが, これは Q。の空間相関を異方的にする. たとえば方向と 方向の相関が抑えられる. これを勾配ベクトル $m$ の相関で言うと, $m$に平行およひ垂直な方向の相関が抑えられることを意味している. 一方 $m$ に対して斜め方向の相関は比較的長い距離まで持続する. このような, 各点での $m$ の方向$|q_{-}qjQ-j|^{2}$$aQ_{ex}$$x$$y$に対して相対的な異方性を定量化するため, 次の相対方向相関関数を定義する.
$\mathrm{o}G_{\mathrm{e}1}(X,\mathrm{Y})=\langle Q_{ij}(r)Q_{\dot{\iota}j}(r+R)\rangle$,$R=(X\cos\theta(r),\mathrm{Y}\sin\theta(r))$

,

(12)

は上記したように $m(r)$ の角度である. この相関関数の $X$ 軸およひ
軸に沿った値はそれそれ, 局所的な勾配ベクトルに対して平行およひ垂直な方向の相関を記述する. 弾性膜と流体膜のそれそれについて, 数値シミュレーションの結果軸上ては, 相から計算した相関関数を図 3 に示す. 弾性膜については $X$ およひ
関距離の約 2\sim 4 倍の距離にあたる領域で相関が負になっている一方, 斜め方向で
ただし$\mathrm{Y}$$\theta(r)$$\mathrm{Y}$は相関はすべての距離にわたって正である. これは上記の弾性相互作用の異方性から予想される傾向と一致する. 対照的に, 流体膜についてはすべての方向で振動的な相関が見られる.

136(b)(a)

図 4: 基板上の膜の座屈パターン. (a) ソフト基板上のミクロ座屈. (b) ’‘-b. 基板からの剥離. いすれも粗大化過程における高さ のスナップショット. グレースケールは (&) では明色が $h>0$ , 暗色が $h<0$ に, (b) では $h=0$ が黒に対応.
$h$5

基板上の弾性膜の座屈前節までの結果から分かるように, 単独の膜の揚合, 平衡状態における座屈パターンの特微的な長さは境界条件のみによって決まる. バターン自体も境界の形状に強く依存するため, 理論と実験の比較の上では単独膜は必すしも理想的な系とは言えない. 一方, 基板によって支持された膜を座屈させた場合, 基板の変形がパターンの成長を止めるか遅くするため, パターンの統計的な特徴が比絞的よく観察されている. また, 支持膜の場合, 座屈は基板と膜の熱膨張係数の差を利用して, 温度変化によって き起こすことができる. 高温で接着した基板と膜を冷却したとき, 基板の収縮率が膜より大きければ, 基板は膜に水平方向の圧力を及ぱすからである (バイメタルの原理) このような実験系の例として (i) ゴムのようなソフト基板上の膜のミクロ座屈 (ii) 金属, ガラスのようなハード基板上の膜の剥離の 2 つが挙けられる. (i)の典型的な例ではエラストマー基板に厚さ数十の金属ナノ薄膜を蒸着したものが用いられる. 膜は基板に接着したまま座屈変形し ,オーダーの幅を持つ皺からなる迷路状パターンが得られている [8]. (ii) の例ではオーダーの厚さの金属蒸$\mathrm{B}1$$|$$\mathrm{n}\mathrm{m}$$\mu \mathrm{m}$$\mathrm{m}\mathrm{m}$

着膜が座屈して基板から剥離し , 剥離した部分は細い網目状につながったバターンを形成する [9].
このような基板上の座屈の最も簡単なモデルとして, 膜と基板との相互作用を垂直ばねで置き換えたものを考える. 基板の変形や剥離の非可逆性は考えないものとして, 前節までの膜のモデルに, 高さに依存するボテンシャル $V(h(r))$ を加える. 具体
的には, ソフト基板の場合, 線形ばねを仮定して$V(h)= \frac{\gamma}{2}h^{2}$,(13)ハード基板の場合, 線形ばねが臨界長 h。で切れると仮定して,$V(h)= \frac{\gamma}{2}(h^{2}-h_{e}^{2})\theta(h_{\epsilon}-h)$(14)

137は階段関数) それそれに対してシはばね定数,ミュレーションから得られたパターンを図 4 に示す、 (i) のミクロ座屈におけるストライプの幅 は, ばねのエネノレギーと曲け弾性のバランスから,なるポテンシャルを用いた ($\gamma$$\theta(x)$-$w\sim(\kappa/\gamma)^{1/4}/\sqrt{|U|}$$w$と見積もられる. ストライプの向きの相関距離は時間とともに成長し,を超えると前節で述べた有効相互作用が支配的となって異方的相関が顕著となる. 一方 (ii)$\sqrt{\gamma/\mu}$の剥離においては核生成型の時間発展が見られた. 初期条件として斑点状の剥離パターンを与えたところ, 大きな斑点だけが線状に伸長し , その再末端に 字型分岐が形成された. この構造は, 前節でのリッジに相当する. その後, このリツジ状ドメインが伸長, 連結して全体としては網目状の剥離パターンとなった. 網目のサイズや網の大さは時間と共に増大するが, その発展則については今後の課題である. また,$\mathrm{Y}$

実験て観察されている蛇行したリッジの成因についても検討の余地が残されている.

6 結論

結果は今後, 大変形およひ不可逆変形を扱う上で基礎になると思われる.
本研究は科学研究費補助金 [ソフトマターにおける方向自由度と粘弾性の協同効果」を受けて行なわれた. 本研究の詳細は論文 [10] およひその他に発表予定である.

参考文献
Treatise on the Mathematical Theory of[1] A. E. H. Love,Edition (Dover, New York, 1944).
$u\mathrm{A}$$\mathrm{E}1\mathrm{a}\mathrm{s}\mathrm{t}\mathrm{i}\mathrm{c}\mathrm{i}\mathrm{t}\mathrm{y}^{n},$$4\mathrm{t}\mathrm{h}$[2] T. A. Witten and H. Li, Europhys. Lett. 23, 51 (1993); A. E. Lobkovsky et ,Science 270, 1482 (1995); E. M. Kramer and T. A. Witten, Phys. Rev. Lett.
Lett. 87,78, 1303 (1997); B. A. Didonna and T. A. Witten, Phys.
(2001).
$a\mathrm{F}.$$\mathrm{R}\epsilon \mathrm{v}.$[3] D. Moldovan and L.
$\mathrm{E}00,4377$$\mathrm{G}\mathrm{o}\mathrm{l}\mathrm{u}\mathrm{b}\mathrm{o}\mathrm{v}\mathrm{i}\acute{\mathrm{c}}$$2\mathrm{o}\mathrm{e}105$, Phys. Rev. Lett. 82, 2884 $(1\Re 9)$ ;Phys. Rev.

(1999).

[4] M. Ben Amar and Y. Pomeau, Proc. R. Soc. Lond. A 453, 729 (1997).

138[5] S. Chaieb and F. Melo, Phys. Rev. Lett. 80, 2354 (1998); E. Cerda and L.
Mahadevan, Phys. Rev. Lett. 80, 2358 (1998); E. Cerda et $al.,$ Nature (London)
401, 46 (1999); E. Cerda and L. Mahadevan, Phys. Rev. Lett. 90, 074302-1(2003).

[6] A. Boudaoud et $al.,$ Nature (London) 407, 718
Boudaoud, Phys. Rev. Lett. 91, 086105-1 (2003).

$(2\alpha 10)$;B. Audoly and A.

[7] L. D. Landau and E. M. Lifshitz, Theory of Elasticity (Pergamon, Oxford,1986).
[8] For areview, aee: G. Giola and M. Ortiz, Adv. Appl. Mech. 33, 119 (1997).

[9] N. Bowden et $al.,$ Nature (London) 393, 146 (1998); N. Bowden et $d.$ , Appl.
Phys. Lett. 75, 2557 (1999); D. B. H. Chua et $d.$ , Appl. Phys. Lett. 76, 721(2000).

[10] N. Uchida, Physica$\mathrm{D},$in press.