計算数理工学論文集 Vol.21 (2021 年 12 月), 論文 No.09-211218

JASCOME

ヘリンボン及びトリプルジャンクションパターンの発達に対する影響因子解析
EFFECTIVE FACTORS ANALYSIS ON EVOLUTION OF HERRINGBONE AND TRIPLE-JUNCTION PATTERNS田中 大地 1)，菊池 正太郎 2)，松原 成志朗 3)，永島壮 4)，奥村大 5)

Daichi TANAKA, Shotaro KIKUCHI, Seishiro MATSUBARA, So NAGASHIMA and Dai OKUMURA
1) 名古屋大学大学院工学研究科（〒464-8603 名古屋市千種区，E-mail: tanaka.daichi@e.mbox.nagoya-u.ac.jp）
2) 名古屋大学大学院工学研究科（〒464-8603 名古屋市千種区，E-mail: kikuchi.shotaro@f.mbox.nagoya-u.ac.jp）
3) 名古屋大学大学院工学研究科（〒464-8603 名古屋市千種区，E-mail: seishiro.matsubara@mae.nagoya-u.ac.jp）
4) 名古屋大学大学院工学研究科（〒464-8603 名古屋市千種区，E-mail: so.nagashima@mae.nagoya-u.ac.jp）
5) 名古屋大学大学院工学研究科（〒464-8603 名古屋市千種区，E-mail: dai.okumura@mae.nagoya-u.ac.jp）
In this study, we analyze effective factors on evolution of herringbone (HB) and triple-junction (TJ)
patterns occurring on a compressed film bonded to a soft substrate, which in turn is bonded to arigid support. Three-dimensional finite element analysis is performed assuming the occurrence ofthe HB and TJ patterns, which are caused by superposition of the hexagonal dimple mode andstripe modes in three symmetric directions. Pattern evolution relatively decreases the total elastic
energy of film/substrate system as the overstress increases, obtaining the quantitative difference
between the HB and TJ patterns by parameterizing film and substrate properties, i.e., the Young’s
modulus ratio (Ef/Es=10~1000) and the thickness ratio (hs/hf =1~100). It is found that the HB
pattern is dominant in most cases but the TJ pattern becomes comparable in a specific case of asmall Young’s modulus ratio and a small overstress. The decrease of hs/hf makes the HB patternsuperior to the TJ pattern.
Key Words: Pattern evolution, Surface instability, Buckling, Wrinkling, Finite element analysis
1. 緒言

表面不安定(1)~(3)は，固体力学の分野において，古典的ではあるが挑戦的な研究課題として近年盛んに研究されており，哺乳類の大脳皮質のしわの形成に代表されるように(4)~(7)，生体システムの発達や形態形成の基礎的理解と結び付くことが期待されている．したがって，軟質基盤上に硬質膜という組み合わせを最も基本的な構成要素とみなし，硬質膜に作用する面内圧縮応力を駆動力として，パターンの発生・発達を研究することが重要となる(8)~(11)．
パターン化の初期段階では，膜表面にはディンプルと呼ばれるくぼみが周期的に生じることが実験的に観察されており(12),(13)，理論的には臨界応力(2)においてリンクルと呼ばれる

Fig. 1 Schematic illustration of surface patterns: (a) herringbone(HB) pattern and (b) triple-junction (TJ) pattern (17).

正弦波形のしわが生じ，このリンクルは異なる方向に任意に生じ得ることから，その重畳発生によってディンプル構造は実験では，ヘリンボンパターンが観察されることが多いた形成される(8),(9),(11)．一方，この臨界応力を基準としてさらに応力が増大するにしたがって，パターンは異なる形態に発達

め，このパターンに着目して研究は進められたが(8),(9),(11)，超する．最終的に大きな過応力状態において，ヘリンボンやラ臨界状態において膜表面にランダムな摂動負荷を与える有ビリンスパターン（ランダムヘリンボン）が観察されており，限要素解析(14)によって，ヘリンボンパターンとカウンターパ中間的な状態でも多様なパターンが生じ得る(12),(13)．

ートとなるトリプルジャンクションパターンの存在が明らかとなった（Fig. 1）．この解析では，基盤と膜のヤング率比2021 年 10 月 6 日受付，2021 年 11 月 10 日受理

が 10 以下の場合を想定しており，この領域においてトリプ− 51 −

ルジャンクションパターンがヘリンボンパターンとエネルギー的に共存可能となり，大脳皮質に生じるしわの形態にはこの共存の機構が重要であることが議論されている(14)．しかし，共存関係に影響する因子の解析は十分ではなく，ランダムな摂動負荷を用いているため，ヘリンボンやトリプルジャンクションパターン自身の発生機構も明らかではなかった．
一方，著者らの研究グループでは，有限要素解析において，逐次的な座屈固有値解析(15)を実施することによって，ディンプルパターンからヘリンボンパターンに発達する過程にて生じる逐次的な分岐点を同定し，どのような分岐モードの発生によって変形パターンが発達するかを明らかにする解析を進めている(16),(17)．結果として，ヘリンボンパターンは第Fig. 2 Initial dimensions of a compressed film bonded to asoft substrate, which in turn is bonded to a rigid support.

一分岐点で生じるディンプルモードと第二分岐点で生じるストライプモードや長方形チェッカーボードモードによって形成されることを示した(17)．さらに，第二分岐点は，ディンプルパターンの対称性に起因して多重分岐となり，分岐モードの選択的な重畳によってトリプルジャンクションパターンも生じ得ることを示した(17)．すなわち，ヘリンボンパターンとトリプルジャンクションパターンは同じ分岐点から派生する兄弟パターンであり，この知見を用いれば，座屈後解析によって系統的にパターン発達の解析を行い，結果を比較することができる．
そこで本研究では，ヘリンボン及びトリプルジャンクションパターンの発達に対する影響因子を有限要素解析によって幅広く明らかにすることを考える．2 章では，ヘリンボン

Fig. 3 Bifurcation modes of (a) hexagonal dimple mode(1)
( HEX) at the first bifurcation and (b)–(d) stripe modes in three(2)symmetric directions ( 0(2) , +60and (2)60 ), respectively, atthe second bifurcation on the first bifurcated path(17).

及びトリプルジャンクションパターンが生じることを前提とする簡易的な解析方法を提案する．3 章では，有限要素モデルについて述べる．解析では，膜と基盤のヤング率比及び膜厚比をパラメータとする．4 章では解析結果を示す．パターン発展過程のエネルギー差を詳細に解析することによっ生する兄弟パターンであることが明らかにされており，Fig.て，大半の場合にヘリンボンパターンが優勢であることと，3 に第一及び第二分岐点で生じる分岐モードをそれぞれ示す．

特定の小さなヤング率比と小さな過応力状態の組み合わせ

図中のコンターは分岐モードの主成分である面外方向成分においてトリプルジャンクションパターンがヘリンボンパに対応している．第一分岐点では，六方配列のディンプルモターンと同等になり得ることを示す．さらには，膜厚比の減(1)ード（ HEX）が生じることを想定しており(11)~(13)，このため少はヘリンボンパターンを優勢にすることも示す．最後に 5L2/L1=31/2 の値を用いている．第二分岐点では，六方ディンプ

章では結言を述べる．

ルモードの対称性に起因して，等価な 3 つの方向にストライ(2)プモード（ 0(2) ， +60， (2)60 ）がそれぞれ生じる．したがっ

2. 解析方法て，第二分岐点では 3 つの分岐モードの組み合わせに依存し本研究で用いる解析方法は，Kikuchi らの研究(17)によってて異なる分岐経路が現れる．ヘリンボンパターンは 0(2) （対得られた知見に基づいており，Fig. 2 に膜基盤構造体の代表(2)
称性を考えると，+60もしくは (2)，トリプルジャ60 でも同じ）体積要素を示す．膜と基盤の層間は完全接着されており，膜(2)ンクションパターンは 0(2)  +60 (2)60（符号の組み合わせにの自由表面（x3=hs+hf）は応力フリー，基盤の底面（x3=0）によらない）が作用することによって，第二分岐点からそれぞは完全拘束を与える．また，側面には周期境界条件を与える．
なお，膜に生じる面内圧縮応力は熱ひずみによって与え，温れパターン発展する．Fig. 4 は対応する分岐モードが重畳されたコンター図であり，発達した後のヘリンボン及びトリプ度上昇を負荷パラメータとする．有限要素解析では，膜厚 hf

ルジャンクションパターン（Fig. 1）の特徴をおおよそ表しで無次元化した寸法を考え，膜厚比 hs/hf の他に L1/hf と L2/hfていることがわかる．

が設定すべき寸法となる．したがって，hs/hf の値が十分に大Fig. 3 と Fig. 4 は特定の寸法や材料定数を設定して座屈固きければ，基盤底面の完全拘束の影響は無視できるようにな有値解析を実施することによって得られた結果であるが，本り，小さくなるにつれてこの影響は現れる．本研究では，膜研究では座屈固有値解析をせずに，ヘリンボンパターンとト厚比の影響をパラメトリックに考える（3 章参照）．

リプルジャンクションパターンの発達する分岐経路を解析

Kikuchi らの解析(17)によって，ヘリンボンパターンとトリすることを考える．このため，ヘリンボンパターンを解析す

プルジャンクションパターン（Fig. 1）は同じ分岐点から派− 52 −(1)
 0(2)
Fig. 4 Initial imperfection consisting of (a) HEX
(1)
(2)
(2)
evolving to herringbone pattern and (b) HEX  0  +60  (2)60
evolving to triple-junction pattern(17).

(1)
 0 (2) and
Fig. 6 Initial imperfection consisting of (a) HEX
(1)
(2)
(2)
(2)
(b) HEX  0  +60  60 (cf. Fig. 4).

Fig. 5 Bifurcation modes consisting of superposition ofsimple sinusoidal waves with amplitude in x3 direction: (a)
(1)
(2)
(2)
HEX, (b) 0 (2) , (c) +60
, and (d) 60
(cf. Fig. 3).

Fig. 7 Finite element meshes of representative volumeelement.

(1)
るときには，初期不整として HEX 0(2) を用い，トリプルジャンクションパターンを解析するときには，同様にして，それぞれ 122,330，115,200 であり，膜は厚さ方向には 10 分(1)(2)
HEX 0(2)  +60 (2)60 を初期不整として用いる．なお，第二分割，基盤は厚さ方向に 15 分割されている．2 章で述べたよう岐点の分岐モードは，第一分岐経路上で生じる六方ディンプ(1)に，六方ディンプルパターンのためには HEX，ヘリンボンパルパターンの影響を受けて，若干ゆがんでいることがわかる(1)
 0 (2) ，トリプルジャンクションパタ
ターンのためには HEX

（Fig. 3）．しかしながら，初期不整として用いる分には，こ(1)(2)(2)ーンのためには HEX 0 (2)  +60 60で構成される初期不整のゆがみを考慮する必要はないと考えられ，本研究では，Fig.

を与える．したがって，膜表面の節点にのみ初期不整を与え，(1)5 に示されるように，六方ディンプルモード（ HEX）及びス

その大きさは，膜厚 hf で無次元化し，面外方向成分の最大値(2)(2)
トライプモード（ 0 (2) ， +60， 60）を x3 方向にのみ振幅をを 0.01 に規格化して与える(16),(18)．

有する正弦波によって生成し，膜表面にのみ小さな大きさで膜には等方的な線膨張係数を=0.01 として与え，温度 T与える．この結果，Fig. 5 に示される簡易モードからヘリンを上昇させることによって，面内圧縮応力を増加させる．線

ボン及びトリプルジャンクションパターンを誘起する初期膨張係数や温度は本質的には意味のない値であるから，パタ不整は Fig. 6 のように作成される．したがって，寸法や材料ーン発展を伴わない主経路上で得られる値を用いて無次元定数を変更しても，座屈固有値解析によって分岐モードを繰化する(8),(11),(17)．すなわち，主経路上では膜の面内に生じるり返し解析する必要性を伴わない．

応力は均質な等二軸圧縮状態となるため，温度に対応した x1及び x2 方向の圧縮応力を  0 (T ) (>0)，さらには第一分岐点で3. 解析モデルの臨界圧縮応力を  0c (Tc ) と表し，以後，T や Tc を用いず，過有限要素解析に用いる要素分割を Fig. 7 に示す．第一分岐

応力状態を説明する上で便利な  0 /  0c の値を用いる（以後，点において六方ディンプルモードの発生を想定するため，この値を過応力と呼ぶ）．同様にして，パターン発展に対す

L2/L1=31/2 となっており，膜厚比には十分に大きな値としてるエネルギー的な比較のため，代表体積要素中に蓄えられるhs/hf =95 を設定する．L1/hf（もしくは L2/hf）の値は後述する

弾性ひずみエネルギーを U(T)と表し，主経路上での値 U0(T)ように同定する．ただし，この値に関わらず，Fig. 7 に示さを用いて無次元化する．したがって，第一分岐点を超えてパれる要素分割を用いる．Abaqus で解析することを考え，三ターン化が進むと（すなわち  0 /  0c >1），U/U0 は減少し，そ次元ソリッド要素 C3D8RH を採用する．節点数と要素数はの程度を比較することによって，発達するパターンの優劣を− 53 −

評価する．
本研究では，ヤング率比と膜厚比をパラメトリックに考えて影響因子を調べる．膜と基盤は，非圧縮性を有する
neo-Hookean 超弾性モデルに従うと仮定し，ポアソン比は 0.5となる．ヤング率をそれぞれ Ef，Es として，ヤング率比 Ef/Es=10, 100, 1000 の場合を解析する．一方，膜厚比は，用意した代表体積要素の底面寄りの節点を完全拘束することによって，hs/hf <95 の場合も解析する(18)．節点の位置に依存するため，hs/hf の値は必ずしもキリの良い数字にはならないが，改めてモデルを作成する必要がないため便利である．
最後に，L1/h（もしくはL2/hf）の決定方法について述べる．
fFig. 8

第一分岐点で生じるリンクルの波長は，理論解(2),(8)に従うこ

Comparison of  0c /Ef as a function of hs/hf.とが知られており(16)~(18)，波長th はth/hf =2(3Es/Ef)1/3，臨界圧縮応力th はth/Ef=(1/3)(3Es/Ef)2/3 と表される．したがって，六方ディンプルモード（Fig. 3(a)，Fig. 5(a)）を構成するリン
クル波形の波長がth/hf となるように L1/hf の値を決定することができる．上述の理論解は，膜と基盤は非圧縮性を有すると仮定し，膜厚比は十分に大きい場合（hs/hf→∞）の結果である．膜厚比が有限の値を取り，小さくなる場合の理論解(10),(18)も存在しており，この場合には，ある閾値を超えると波長は短くなり，臨界圧縮応力は増加する傾向を有する．本
研究では，影響因子の解析を主目的とし，解析結果の比較を容易にするために，膜厚比が波長に及ぼす影響を考慮せず，

 0c   th として過応力  0 /  0c を評価する．
4. 解析結果(1)
Fig. 8 は，六方ディンプルモード HEXを初期不整に用いるときの，臨界圧縮応力のヤング率比及び膜厚比に対する依存性を示している．マーカー付きの実線は有限要素解析によって得られた値，鎖線は hs/hf→∞の場合の理論値である．なお，有限要素解析では U/U0 が 1 から 0.1%減少したところの圧縮応力の臨界値とみなした．この図が示すように，膜厚比が大きい場合（hs/hf =95）には，解析値は理論値と良い一致を示しており，2・3 章で述べた解析手法とモデル化がうまく機能

Fig. 9 Comparison of U/U0 as a function of  0 /  0c for (a)Ef/Es=10 and (b) Ef/Es=1000.

している．膜厚比が小さくなると，徐々に基盤底面から拘束効果を受けるようになり，プラトー領域を経て，ヤング率比に依存して臨界圧縮応力は増加傾向を示すようになる．臨界なわち，ヤング率比が高く，膜厚比が小さい条件下ではヘリ圧縮応力が変化すると，リンクルの波長も変化すると考えらンボンパターンが明らかに優位となることがわかった．

れるため(10),(18)，本解析手法では図中にプロットしており，Fig. 10 と Fig. 11 には，  0 /  0c =2.5 における面外方向の変臨界圧縮応力が増加傾向を示す前までの膜厚比を解析範囲

位のコンター図を示しており，ヤング率比及び膜厚比がヘリとみなし，ヘリンボンパターンとトリプルジャンクションパンボン及びトリプルジャンクションパターンの発展に及ぼターンの発現優位性について比較する．

す影響を示している．これらの図が示すように，ヘリンボンはじめに，ヘリンボン（HB）及びトリプルジャンクショパターンの発達には，ヤング率比や膜厚比の影響がほとんどン（TJ）パターンの解析によって得られるエネルギー変化と現れないことがわかる（Fig. 10）．したがって，パターン発過応力の関係を Fig. 9 に示す．この図は，Ef/Es =10, 1000 の達によるエネルギー減少にもこれらの変化は陽に現れない場合であり，膜厚比が十分に大きく，基盤底面からの拘束効と考えられる．一方，トリプルジャンクションパターンの発果が無視できる場合には，ヘリンボンパターンとトリプルジ達には，これらの影響は顕著であることがわかった（Fig. 11）．

ャンクションパターンの間のエネルギー差はほとんど優劣

すなわち，ヤング率比が大きくなるとトリプルジャンクショのないようにみえる．一方，膜厚比が小さくなると（hs/hf

ン部の変形が幅広となり，さらに膜厚比が小さくなると，そ=28.5），ヤング率比が高い場合（Ef/Es =1000），トリプルジャの中心部分の面外方向の変形が抑制される．したがって，パンクションパターンのエネルギー減少は鈍くなる．ヘリンボ

ターン発達によるエネルギー解放が抑制されることとなり，ンパターンのエネルギー減少にその傾向は現れておらず，す− 54 −

Fig. 10 Deformed configurations at  0 /  0c =2.5 in the caseof herringbone pattern with (a) {Ef/Es, hs/hf}={10, 95}, (b){Ef/Es, hs/hf}={10, 28.5}, (c) {Ef/Es, hs/hf}={1000, 95}, and (d){Ef/Es, hs/hf}={1000, 28.5}.

Fig. 12
Energy difference, (UTJUHB)/U0, between
herringbone and triple-junction patterns as a function of hs/hf
and  0 /  0c ; (a) Ef/Es =10, (b) Ef/Es =100, and (c) Ef/Es=1000.
12 に示す．なお，UHB と UTJ はヘリンボンパターンとトリプ
ルジャンクションパターンが発達するときの代表体積要素中に蓄えられる弾性ひずみエネルギーを示している．したがって，(UTJUHB)/U0 の値が正の場合にはヘリンボンパターン，負の場合にはトリプルジャンクションパターンが優勢となり，その差が大きくなるほど定量的にもこの傾向が強いと理Fig. 11 Deformed configurations at  0 /  0c =2.5 in the caseof triple-junction pattern with (a) {Ef/Es, hs/hf}={10, 95}, (b){Ef/Es, hs/hf}={10, 28.5}, (c) {Ef/Es, hs/hf}={1000, 95}, and (d){Ef/Es, hs/hf}={1000, 28.5}.

解できる．
これらの図を見ると，膜厚比が十分に大きく，基盤底面の拘束効果が無視できる場合には，過応力が大きくなると，ヤング率比に関わらず，エネルギー差は 0.005 あたりに分布しこのような機構によって，ヘリンボンパターンが優勢となり，ている．一方，ヤング率比が小さい場合（Ef/Es=10）には，トリプルジャンクションパターンが劣勢になると推察でき過応力レベルが小さい初期段階（  0 /  0c =1.5）において，エる．

ネルギー差は負の値を取り，ヘリンボンパターンではなくトより詳細な定量的検討のため，ヘリンボンパターンとトリリプルジャンクションパターンが陽に優勢となることがわ

プルジャンクションパターンのエネルギー差の変化を Fig.

かった．すなわち，この過応力レベルの周辺が，文献(14)で議論されているヘリンボンパターンとトリプルジャンクシ− 55 −

ョンパターンの共存領域に対応すると考えられる．このときのエネルギー差の絶対値は 0.001 程度であるから，過応力が増加してエネルギー差が正の値として 0.005 を取る場合というのは，ヘリンボンパターンが優勢であると類推でき，文献(14)の結果と合わせて整合的に解釈することができる．
本研究では膜厚比の影響も解析しており，膜厚比が小さくなると，エネルギー差は正の方向に拡がり，すなわち，先に述べたようにヘリンボンパターンが優勢となる．エネルギー差は，ヤング率比によって変化し，図中では 0.010～0.020 まで過応力の増加に伴って増加する．すなわち，基盤底部の拘束効果はヘリンボンパターンを優勢にすることがわかった．
5. 結言

本研究では，軟質基盤上の硬質膜に生じるヘリンボンパターンとトリプルジャンクションパターンの発達に及ぼす影響因子の有限要素解析を行った．座屈固有値解析を実施しな
い簡易的な解析方法を用いたが，どちらのパターンについてもその発達過程を解析することができた．詳細なエネルギー差の比較によって，大半の場合においてヘリンボンパターンが優位であり，特定の小さなヤング率比と小さな過応力状態の組み合わせにおいてトリプルジャンクションパターンが同等もしくは優勢であることがわかった．この知見は文献(14)の結果と整合する．また，膜厚比が小さくなると，基盤底面の拘束効果が生じるようになり，この効果はヘリンボンパターンを優勢にすることがわかった．この原因として，トリプルジャンクションパターンの変形はこの拘束効果によって抑制されるのに対して，ヘリンボンパターンにはこの抑制機構が現れないためであると考えられる．
最後に，今後の課題について述べる．本研究では簡単のため，ポアソン比が 0.5（非圧縮）の場合を想定したが，圧縮性がパターン発展に及ぼす影響を調べることは重要である．
また，エネルギー差を評価してパターン発展の優位性を議論したが，このエネルギー差がパターン発展の安定性に及ぼす影響については異なる議論が必要かもしれない．さらに，パターンが発達すると，自己接触を伴うクリースが発生する可能性があり，クリースとリンクルの相互作用解析は新しく必要な解析技術となり得る．このような解析技術の発達と解析結果の蓄積によって，生体システムに生じる複雑なパターンの固体力学的理解が進展すると考えられる．
謝辞
本研究は JSPS 科研費 JP19H00739 の助成を受けて行われた．

ここに記して謝意を表する．
参考文献(1) M.A. Biot: Mechanics of incremental deformations: theory ofelasticity and viscoelasticity of initially stressed solids andfluids,including
thermodynamic
foundations
andapplications to finite strain, 1965, Wiley.
(2) H.G. Allen: In: Analysis and design of sandwich panels,Pergamon Press, 1969, New York.

− 56 −(3) A.N. Gent, I.S. Cho: Surface instabilities in compressed orbent rubber blocks, Rubber Chemistry and Technology,
72(1999), pp.253–262.
(4) D.P. Richman, R.M. Stewart, J.W. Hutchinson, V.S. CavinessJr.: Mechanical model of brain convolutional development,
Science, 189(1975), pp.18–21.
(5) S. Budday, P. Steinmann, E. Kuhl: The role of mechanics
during brain development, Journal of the Mechanics andPhysics of Solids, 72(2014), pp.75–92.
(6) T. Tallinen, J.Y. Chung, F. Rousseau, N. Girard, J. Lefèvre, L.
Mahadevan: On the growth and form of cortical convolutions,
Nature Physics, 12(2016), pp.588–593.
(7) L.da.C. Campos, R. Hornung, G. Gompper, J. Elgeti, S.
Caspers: The role of thickness inhomogeneities inhierarchical cortical folding, Neuroimage, 231(2021),
No.117779.
(8) X. Chen, J.W. Hutchinson: Herringbone buckling patterns ofcompressed thin films on compliant substrates, Journal ofApplied Mechanics, 71(2004), pp.597–603.
(9) B. Audoly, A. Boudaoud: Buckling of a thin film bound to acompliant substrate—Part I: formulation, linear stability ofcylindrical patterns, secondary bifurcations, Journal of theMechanics and Physics of Solids, 56(2008), pp.2401–2421.
(10) Z.Y. Huang, W. Hong, Z. Suo: Nonlinear analyses of
wrinkles in a film bonded to a compliant substrate, Journal of
the Mechanics and Physics of Solids, 53(2005),pp.2101–2118.
(11) S. Cai, D. Breid, A.J. Crosby, Z. Suo, J.W. Hutchinson:
Periodic patterns and energy states of buckled films oncompliant substrates, Journal of the Mechanics and Physicsof Solids, 59(2011), pp.1094–1114.
(12) D. Breid, A.J. Crosby: Effect of stress state on wrinkle
morphology, Soft Matter, 7(2011), pp.4490–4496.
(13) M. Guvendiren, S. Yang, J.A. Burdick: Swelling-induced
surface patterns in hydrogels with gradient crosslinking
density, Advanced Functional Materials, 19(2009),
pp.3038–3045.
(14) T. Tallinen, J.S. Biggins: Mechanics of invagination andfolding: hybridized instabilities when one soft tissue grows
on another, Physical Review E, 92(2015), No.022720.
(15) D. Okumura, J. Sugiura, H. Tanaka, Y. Shibutani: Bucklingand postbuckling of etching-induced wiggling in a bilayer
structure with intrinsic compressive stress, International
Journal of Mechanical Sciences, 141(2018), pp.78–88.
(16) H. Miyoshi, S. Matsubara, D. Okumura: Bifurcation anddeformation during the evolution of periodic patterns on a gel
film bonded to a soft substrate, Journal of the Mechanics andPhysics of Solids, 148(2021), No.104272.
(17) S. Kikuchi, S. Matsubara, S. Nagashima, D. Okumura:
Diversity of the bifurcations and deformations on films
bonded to soft substrates: robustness of the herringbone
pattern and its cognate patterns, (submitted).
(18) H. Miyoshi, D. Okumura: Effect of soft substrate thickness
on wrinkle analysis of a gel film, Transactions of the Japan
Society for Computational Methods in Engineering, 19(2019)，
pp.31–36 (in Japanese)．