粘着剤の設計支援に向けた
Euler型有限要素法に関する研究( Eulerian finite element methodfor design of pressure sensitive adhesive

学位取得年月(平成22年9月)
前田 和久

目次
第1章 緒言1.1背景と目的
1.2 本研究の概要

第2章 粘着剤の力学特性
2.1粘着テープ.
2.2 粘着剤の変形.
2.2.1 ゴム弾性‥

82.2.2 粘弾性

112.2.3 ピール試験での変形

15

2.3 粘着剤の解析

第3章

16

2.3.1粘着特性の解析

172.3.2 材料モデル

18

2.3.3 シミュレーションに関する既存の研究

19

Euler型有限要素法

20

3.1保存則20

3.1.1質量保存の法則

20

3.1.2 運動量保存の法則

223.2 0peratorSplit法.

23

3.3 非移流ステップ.

24

3.3.1 平衡方程式の離散化

243.3.2 動的陽解法.

25

3.3.3

アイソパラメトリック要素‥

3.4 移流ステップ
3.4.1 密度関数‥

. 27

29
. 29

3.4.2 速度

30

3.5 VOF法による界面捕捉

30

3.5.1 VOF法

30

3.5.2 2次精度風上差分法‥. 31

3.6 時間増分制御

36

第4章 粘着剤の構成方程式
4.1微小変形理論における等方性弾性体の構成方程式
4.2 一般化Maxwellモデル

39

4.2.1 ‑般化Maxwellモデルの定式化

40

4.2.2 物質客観性の原理‥。 42

4.3 一般化Maxwellモデルの修正

48

第5章

4.3.1弾性係数‥. 484.3.2 速度形定式化

51

3次元PI」IC法による界面捕捉

54

5.1 3次元PLIC法の概要

55

5.2 境界面の再構成

56

5.2.1境界面セルの判定‥. 56

5.2.2 線形近似された境界面の方程式

565.2.3 境界面の法線ベクトルの計算

57

5.2.4 距離定数の計算.

m15.2.5 境界面の場合分け.

69

5.3 境界面の移流.

735.3.1移流方程式の離散化.

73

5.3.2 体積フラックスの計算

74

第6章 境界条件の付与と空間離散化

75

6.1境界条件の付与.

。 75

6.1.1有限要素法による空間離散化

77

6.1.2 選択型低減積分

78

6.1.3 力学的境界条件の処理

81ii

6.1.4 幾何学的境界条件の処理

83

第7章 解析結果

85

7.1一般化Maxwellモデルの検証

85

7.1.1応力緩和解析

85

7.1.2 ‑事由引張解析

887.1.3 日垂変形解析.

89

7.2 修正した一般化Maxwellモデルの検証

91

7.2.1材料パラメータの取得

92

7.2.2 ‑軸引張試験

937.3 移流計算テスト

96

7.3.1

2D Zalesak's disk rotation

96

7.3.2

2D Shearing月Iow

97

7.3.3

3D Stan ford bunny rotation

7.3.4 計算時間‥. 105

7.4 固体‑固体達成問題

1087.4.1 2個の超弾性体の自由落下問題‥
7.5 固体境界面の‑車由引張問題‥

. 108
. 112

7.5.1力学的境界条件による‑軸引張112

7.5.2 幾何学的境界条件による‑軸引張115

第8章 結言117

付録A非圧縮性固体及び流体の構成方程式

121

A.1超弾性体‥. 121

A.I.I 非圧縮性neo‑Hooke体

121

A.1.2

123

左Cauchy‑Green変形テンソルの評価

A.2 非圧縮性Newton流体

123

謝辞

125

参考文献

126

闇

第1章 緒言
1.1 背景と目的
粘着テープは紙,布,プラスチック,時には金属からなる基材フイルムに,その片面もしくは両面に粘着剤(Pressure Sensitive Adhesive二 PSA)を塗布したもので,材料(被着体)に貼り付けられたり,被着体同士を貼り合わせたりする目的で使用される.つまり,粘着テープの基本となる機能は,貼りつけたり物体同士を貼り合せるというものである.機能としては単純であるが,単純であるがゆえに適応範囲は広い.また,この単純な機能にさらに他の機能を付け加えることで,今までは使用しなかった用途に粘着テープが用いられるようになる.被着体に傷をつけないという機能を付け加えると,表面保護テープという製品となる.家電製品を購入した時に液晶画面に貼り付けられているテープがこれである.一定速度で物質を放出するという機能を付け加えると,治療用テープ製剤となる.これは,例えば職息疾患のある患者の体に貼り付けることで,粘着剤中に含浸させている薬剤を一定速度で体内へ放出し,錠剤等の薬では難しい長時間の安定した薬の効用を実現している.このように粘着テープは基本となる機能は単純であるが,そこに様々な機能を付与することでその適応範囲を広げ,様々な産業分野で用いられるに至っており,またその用途は今でも拡大している.
一方,製品の進歩は目覚しく,粘着テープに要求される機能もますます高くなっている.例えば,近年の電子機器などは軽薄短小と言われるように軽く,薄く,短く,小さくと,より小型化,高密度化が進んでいる.同じ機能を要求されている場
合でも粘着テープとしてはより薄く軽いものが要求される.そして,要求される機能を実現するために,粘着テープの設計はより複雑なものとなってきている.しかし,今日多くの場合,粘着剤の設計は実際に粘着剤を作製し特性を評価することで行われており,その設計には試行錯誤が多い.このようななか,粘着テープの設計に対して数値計算による設計支援が強く要望されるようになっている.

近年,計算機性能の飛躍的向上と計算力学の進歩によって,様々な力学現象の数値解析が可能となっている.この進歩を受けて工業製品の設計に対しても,数値解
析による設計支援は様々な分野で用いられている.自動車や電子機器の設計にはシミュレーションを用いることが日常的に行われ,実際に製品を試作する回数は少なくなってきている.この技術を粘着テープに適応することが可能となれば,粘着テープの設計速度が飛躍的に向上するだけでなく,いまだ判らないことが多い粘着理論に対してもそのメカニズムの解明に貢献でき,今までにない高機能を有する粘着テープの設計が可能となることが期待される.
しかし,粘着テープの挙動を数値解析することは容易ではない.その理由は2つある. 1つは粘着テープの粘着性を発現させる粘着剤は,変形挙動が非常に複雑なことである.もう1つは,粘着剤は小さな力で容易に大変形することである.
粘着剤の応力‑ひずみ関係は単純な弾性の関係にはない.また,金属などの解析においてよく研究されている弾塑性の関係にもない.超弾性に似た非線形挙動を示す.さらに,粘着剤は粘弾性挙動を示す.つまり,非線形弾性の性質以外に応カー
ひずみ関係が時間に依存する粘性の性質も示す.このように変形挙動が非常に複雑であるため,その変形挙動を記述できる適切な材料モデルはない.
また,粘着剤は被着体から剥離されるときに糸曳きと呼ばれる複雑な変形をすることが知られている国14].図1.1に糸曳き現象の模式図を示す.図は上面に被着体があり,下の反っている部分が基材である.その間に粘着剤の変形している様子が描かれている.このように,粘着テープが剥離されるとき,粘着剤の変形は非常に複雑となる.

¥>

Fiff. 1.1: Fingering of PSA

また,粘着剤は容易に大変形するため数値解析が難しくなっている.数値解析手法は,解析格子の制御方法によってLagrange型解法とEuler型解法に分類される.
Lagrange型解法では,図1.2で示すように,物質の変形に追従した観測点において変形を記述する.一方, Euler型解法では図1.3で示すように,空間に固定された解析格子を用いる解法であり,解析格子を超えて物質が変形する.

IA
漢.
+i : I
i>:.
三
t
L.I
LL

Fig. 1.3: Eulerian method

Fig. 1.2: Lagrangian method

通常,国体の数値解析には,経路依存性の問題や固体境界面の明確化という要請からLagrange型解法を用いる.しかし, Lagrange型解法での有限要素法では図1・2のように解析格子が物体に追従して変形するため,極めて大きな変形を取り扱う場
合は計算格子の健全性が損なわれることから解析精度を欠き,さらには計算不能となる.一方,流体解析で主に用いられるEuler型解法での有限要素法では図1.3のよ3

うに解析格子が空間に固定されている.その結果,任意の大変形問題やそれに伴った亀裂や破壊などの現象をLagrange型解法よりも効率的に解明することができる

剛6]閏.
粘着剤の弾性率は室温で105Pa程度であり金属などに比べると極めて低い.もし
粘着斉田こ力が加えられると,それが小さなガであっても粘着斉桐こは極めて大きな変形が生じる.この柔らかいという性質があるからこそ粘着剤は粘着性を発現するのであるが,数値解析を試みると結果として, Lagrange型解法を用いる方法では計算の破綻を起こし解析が困難となる.
本研究では糸曳きに示されるような粘着剤の複雑な変形挙動について数値解析することを最終目標としているが,その第一段階の目標として,粘着剤の大変形解析を可能とすることとしている.この目的を達成するため,粘着剤の大変形解析に関しては,計算格子の破綻が生じないEuler型解法を適応させることで可能とさせる.
また,粘着剤の複雑な変形挙動を記述することに関しては,非線形・大変形を記述できる新しい材料モデルを提唱することとする.また,実質的に粘着剤の変形挙動を解析するためには,物質界面へ境界条件を付与することが必要である.従来のEuler
型解法では困難な物質界面への境界条件の付与について,高精度な界面捕捉法を導入して物質界面を定義し,この定義された物質界面に境界条件を付与することを可能とする.

1.2 本研究の概要
2章では,本研究での取り扱いの対象である粘着剤について説明する.特に粘着剤の変形挙動について,その物性の特徴である非線形,大変形挙動及び応力の時間依存性について,説明する.また,粘着剤の解析に関する今までの検討についても説明する.
3章では,本研究で取り扱うEuler型有限要素法について説明する.本研究で用いるEuler型有限要素法ではアルゴリズムの簡便さから,移流の解法としてOperator
Split法を用いている.また,界面捕捉法としてVOF法を説明する.本研究で用いるVOF法では2次精度風上差分法であるMUSCL法を用いている.

44章では,粘着剤の変形解析を行う上で極めて重要な材料モデルについて説明する.粘着剤の特徴的な性質の1つである粘弾性は一般化Maxwellモデルを用いて記述する.このモデルの定式化においては三次元に拡張して行う.もう1つの性質である非線形,大変形を示すゴム弾性については,一般化Maxwellモデルを修正した
新しい材料モデルを提案する.提案したモデMlついても三次元での定式化を行う.
5章では, Euler型解法で問題である界面捕捉法について説明する.界面捕捉についてはPLIC(Piecewise Linear Interface Calculation)法[8]と呼ばれる手法を用いる.
この手法は流体解析において適用が増えている[9上岡.本研究では3次元PLIC法の構築を行っている.

6章では,境界条件の付与方法について説明する.境界条件は平衡方程式の空間離散化の際に定式に取り入れる.この際,平衡方程式の空間離散化には,有限要素法を用いている.有限要素法を用いることにより,力学的境界条件及び幾何学的境
界条件を取り扱いやすくなる.また,本研究では幾何学的境界条件の定式化にペナ
ルティ法を用いている.
7章では,本手法による解析結果を説明する.まず,粘着剤の変形を記述する材料モデルとして一般化Maxwellモデル及び修正した一般化Maxwellモデルの妥当性を検証する.次いで,移流計算テストの結果を説明し, 3次元PLIC法による界面の捕捉精度を検証する.さらに,従来のEuler型固体達成手法では解析が困難であった,接触を伴う固体一団体達成問題の解析結果を説明する.最後に,固体境界面の‑軸引張問題の解析結果を説明する.これら̀により,本研究にて構築した手法が本研究の目的である粘着剤の大変形解析を可能とすることを示す.
8章では,今までの研究成果をまとめ,粘着剤の大変形解析が可能となったことを示して結言とする.

*

第2章 粘着剤の力学特性
本章では,本研究での取り扱いの対象である粘着剤について説明する.特に粘着剤の変形挙動について,その特徴である非線形,大変形挙動及び応力の時間依存性について説明する.また,粘着剤の解析に関する今までの検討についても説明する.

2.1 粘着テープ
粘着テープは紙,布,プラスチック,時には金属からなる基材フイルムに,その片面もしくは両面に粘着剤を塗布したもので,被着体に貼り付けられたり,被着体同士を貼り合わせたりする目的で使用される.
また,粘着テープの機能は被着体を貼り合わせるだけにとどまらず,保護機能,逮
湿機能,光透過機能,防水・機密機能,導電機能,絶縁機能など様々な機能をさらに付与することで,様々な用途に用いられている.
粘着テープにこれら機能を付与することで使用される用途例を簡単に示すと,セロハンテープやクラフト粘着テープなどの一般事務用テープ,マスキング用テープ,
自動車用塗膜保護フイルムなどの工業用テープ,救急粋創膏,けい皮吸収型テープなどの医療用テープ,熟剥離シート,半導体ウエハ保護・固定用テープなどのエレクトロニクス用テープ,液晶用偏光板などの光学用テープなど多種多様な分野で用いられており,今日なくてはならない製品となっている.
粘着テープの被着体同士を貼り合わせるという機能についてはよく接着剤と対比
古れる.ここで,粘着テープの粘着性を発現させる部材である粘着剤と接着剤について考えてみる.それらが持つ同じ"くっ付く"という機能において異なる点は,くっ付くという機能が発現するまでの時間が挙げられる.
接着剤の場合,機能が発現するには最初液状であった接着剤が硬化反応などを伴い態が液体から固体へと変化することが必要となる.まだ十分硬化していない接着剤で接合された物質は引き剥がそうとすると,非常に弱い力で引き剥がすことがで6

起電口塾瀞

mMmmmsm細管園

クリーン̲ 、̲̲?≠イ1/ j;轟重野 !儀醐

ozsil

問項目間

、屠偉世 子
1

1.、

看■鯵夢野

＼1‑」志一ii〆‑

/

ーI

L

▼.

有権EL LED
こ

r LL

一

一

Fig. 2.1: Applications of PSA

きる.つまり、くっ付くという機能の発現までには時間がかかることになる.一方
粘着剤では,貼った瞬間からくっ付く力が発現される.つまり,機能発現のために時間がかからない(Fig.2.2参照).その反面,接着強さは概して接着剤の方が高い.
しかし,そうであっても時間がかからないという使用上の簡便性や,時には接着強さが弱いということから簡単に剥がして再度貼り合せることができるという再剥離性という性質から粘着剤の使用範囲は広いものとなっている.
接

EV3

強
さ
rf‑^v

m.

付
く

力
ヽー

接着時間(くっ付くまでの時間)
Fig. 2.2: Difference between PSA and adhesive

接着剤は英語で'adhesive"と呼ばれるのに対して,粘着剤は"pressure sensitiveadhesive",日本語に直すと感圧性接着剤と呼ばれる.これは,粘着剤に弱い圧力を7かけるだけで瞬時に被着体とくっ付くという機能が発現されるためである.この粘着剤が弱い圧力でかつ瞬時にくっ付くという機能を発現できる要因は,粘着剤が極めて柔らかいことと粘着剤が単純な弾性という性質ではなく粘性を併せもった粘弾性体であることに起因する16].

2.2 粘着剤の変形
本節では,粘着剤が粘着性を発現する要因である,柔らかく大変形すること,及び弾性と粘性の両方の性質を併せ持つ粘弾性体であることの説明をする.

2.2.1

ゴム弾性

粘着剤の物性を評価する時,しばしば応カーひずみ測定(SS測定)を行う.粘着剤のSS測定は以下のようにして行う.断面積が約2mm?である円柱状の粘着剤を
長さが10mmになるよう引張り試験機に設置する.次いで,一定速度で粘着剤を引き伸ばし,そのときの応カーひずみ関係を測定する. Fig.2.3に一般的な粘着剤のSS曲線を示す.

星
め

誰
a
め
°3

.6
E
CI
l

0%

100%

200%

300%

400%

500%

600%

700%

800%

900%

1000%

nominal strain

Fig. 2.3: SS behavior of PSA

図の横軸は公称ひずみであり,縦軸は公称応力である.粘着剤の弾性率の求め方には様々な方法がある.動的粘弾性試験から求める方法や, ‑軸引張り試験から求8

める方法がある. ‑軸引張り試験から求める場合は, Fig.2.3で示されるSS曲線の初期の立ち上がりにおける接線の傾きから弾性率を求める.
一般的な粘着剤の弾性率は室温で105Pa程度であり金属などの固体材料と比べると極めて低い.そのため,粘着剤に力が加えられると,それが小さな力であっても粘着斉耶こは極めて大きな変形が生じる. Fig.2.3からも,粘着剤は容易に数百%の大変形を生じることが判る.また,種類によっては数千%の大変形を生じることも>‑‑*この粘着剤のSS曲線に見られる非線形・大変形挙動はゴム弾性と呼ばれる性質に起因する.ここで,ゴム弾性の発現機構を説明する.
あらゆる物質は原子や分子の集合体として成り立っている.金属の場合は金属結合によって,また有機化合物の場合は主に共有結合によって原子・分子同士が結びついている.物質を引っ張る(あるいは変形させる)ことは,この原子・分子同士の結合を引き離す作業に他ならない.つまり,分子たちは自分たちにとって安定な状態になるようにお互い結びついて物質を形成している.そこに無理やり引き離そ
うとする力が働くと,それぞれの分子間にはたらく結合力がその抵抗力となり力を発生する.

・・・車●日日●=‥●* ・*‑.・・・.I.・J
A洗er elongated[liiiliulj

Fig. 2.4: Mechanism of enthalpy elasticity

金属などの国体物質の場合,その弾性の発現機構は金属結合や共有結合の結合距
離や結合角度が変化することに起因する(Fig.2.4参照).これは弾性の中でも結晶弾性またはエンタルピー弾性と呼ばれる.金属結合や共有結合の結合距離を伸ばしたり,結合角度を変化させるには大きな力が必要となる.つまり,結晶弾性はその発現のために大きな力が必要となる.また,金属結合や共有結合の変化は大きくはなりえないため,金属などの物質の変形では初期の長さの何倍もの変形は生じない.
9

実際,原子・分子間の結合の歪みによる変形は、初期の長さのわずか1 %にも満たい‥

これに対して,物質の初期の長さに対して何倍もの変形を生じる,粘着剤のSSで見られるような大変形は結晶弾性では引き起こらない.粘着剤のような物質は金属などと異なり高分子と呼ばれる物質から成り立っている.金属は個々の原子が構成単位となり,それらが規則正しく並んでお互いが結合することで形成されているため,それぞれの原子はほとんど自由に動くことができない。これに対し,ゴムのような高分子物質はひとつひとつの分子がとても長く,それらが寄り集まることで作られている。これは,長いひもが適当に絡み合っているような状態のため,規則性
I

は金属に比べ非常に小さく,また互いの結びつきも小さいため,各分子がある程度自由に動き回ることができる.
高分子の存在状態は統計的に考えることができ, Fig.2.5で示すように団子状に折れ曲がった状態になっている.つまり,この状態が高分子にとって最もエネルギー的に安定である.

Fig. 2.5: State of polymer

粘着剤の伸びは,粘着剤を構成するこれらまるまった高分子鎖を一定の方向に引き伸ばすことによって引き起こされる.金属のときのような小さな変化ではなく,粘着剤の場合は,丸まったひもを引き伸ばすという大きな変化のため,初期の長さの何倍もの長さまで伸びることができるようになる(Fig.2.6参照).
最も安定な状態から引き伸ばされると,エントロピーが低下しエネルギーが不安定になる.よって,エントロピー増大則に従って高分子鎖を縮まそうとする力が発生する.この力による弾性のことをゴム弾性またはエントロピー変化によって起こることからエントロピー弾性とも呼ばれる囲.

10

!L‑‑iST‑‑1L

寡一一‑

Remove road

Å鮎r elongatedFig. 2.6: State of elongating polymer

2.2.2 粘弾性
次に粘着剤のもう1つの性質である粘弾性について説明する.
粘着剤は弾性と粘性の性質を併せ持つ粘弾性体である.粘弾性の評価はしばしば
動的粘弾性試験を用いて行われる.動的粘弾性試験とは粘着剤に周期的なひずみを与え(通常正弦波を与える) ,その時の応力の応答を評価するものである.今,式(2.1)で表されるひずみを与えることを考える.
7 ‑ 7osinwf(2.1)

ここで7,7。,のまそれぞれ,ひずみ,ひずみ振幅,周波数である.完全な弾性体であればひずみと応力の関係において位相のずれはない(式(2.2)).また、完全な粘性体であればひずみと応力の関係において位相がtt/2ずれる(式(2.3)).粘弾性では応力とひずみの間の関係は式(2.4)のようになり,位相が∂だけずれる[咽.ここでa,aoはそれぞれ,応力,応力振幅である.
a ‑ cosintot

a ‑ ctqsin (^t ‑ 7172)
a‑ ctqsin(ut‑S)

また,動的粘弾性評価における応力は式(2.5)で示すように複素平面で記述することができる.
a ‑ ao (cosu)t + ismut)(2.5)

このときの弾性率Gサ(複素数で表されるため複素弾性率と呼ぶ)は貯蔵弾性率
G′(W)と損失弾性率G′′u)を用いて式(2.6)のように記述できる.
G*{J)

‑

G′(u)

Ill

+%G′′

u)

(2.(

粘弾性を記述する最も単純なモデルとして, Hooke弾性を示すスプリングとNewton
粘性を示すダッシュポットが直列につながったMaxwellモデルがある.ここで,車ま

よ
(7

n

毒

Fig. 2.7: Maxwell model

粘性係数である.また一次元のモデルを考えているためせん断と引っ張りを考慮せず弾性率をGとしている.
このモデルを用いると貯蔵弾性率G′ u> はその名の通りエネルギーが貯蔵される弾性部であり, Maxwellモデルでいえばスプリング部のGの変化に対応する.一方,
損失弾性率G′′ co はその名の通りエネルギー損失に基づく弾性率であり, Maxwell
モデルでいえばダッシュポット部の77の変化に対応する. Maxwellモデルを考えた
場合 G'(u)とG′′ coは具体的には以下の式で表される[18].
9

G′'u)

‑G

G〝U)

‑G

y

1 +uo2T2

(2.7)

LJ ′「

2^‑2
1+U2T

(2.8)

ここでTは緩和時間と呼ばれV/Gで表される.また,このとき位相のずれ6は式(2.9)で表せる.
tan∂ = 二(2.9)

一般的な粘着剤の動的粘弾性の測定結果をFig.2.8に示す.この測定での周波数は1Hzである.図から粘着剤のtan∂は0ではなく、正の有限の値を示しており,粘弾性を示していることが判る.また, Dahlquistはタックと呼ばれる粘着剤の"ベトベト感"は弾性率が106Pa以下の時に発現されると述べている[19]. Fig.2.8からも粘着剤のG/はタックが発現できるくらい十分小さいことが判る.

12

1.4
s
蝣2

7
1

互6
0.8

召5

‑
00

.̀ロ

星4

0.6

哲3

0.4

2

0.2

0
蝣100

‑50

50

100

150

200

250

Temp. 【℃】

Fig. 2.8: Dynamic mechanical analysis of PSA

粘着剤の粘弾性は動的粘弾性試験だけではなく,応力緩和試験においても見ることができる. Fig.2.9は粘着剤の応力緩和試験の結果を示す.以下に,応力緩和試験の測定方法を説明する.粘着剤のSS測定と同様に,断面積が約2mm?である円柱状の粘着剤を長さが10mmになるよう引張り試験機に設置する.まず,所定の初期ひずみを粘着剤に与える.その後,ひずみ一定にした状態での応力と時間の関係を評価する. Fig.2.9に一般的な粘着剤の応力緩和試験の結果を示す.

200

400

600

800

1000

1200

1400

Time [sec]Fig. 2.9: Stress relaxation behavior of PSA

これから,粘着剤の応力は時間とともに変化する,つまり時間依存性があることが判る.この応力緩和現象は粘弾性が示す典型的な性質の例である.
13

なぜ粘着剤が粘弾性を示すかについて,応力緩和性に基づいて説明する.
高分子材料は粘弾性を示すものが多い.中でも,粘着剤は室温において著しい粘弾性を示す.粘着剤を構成する1つの高分子は他の高分子と共有結合によって結合(架橋と言う)されていないものが多い(一部架橋を施しているが,粘着剤を構成する全ての高分子を架橋すると,粘弾性挙動が失われ粘着特性が低下するため通常は行わなv¥).このように共有結合によって結びつけられていない材料を引き伸ばす際に生じる抵抗力は, Van der Waals力や水素結合などに起因する分子同士の親和力以外に高分子の絡まりあいにより変形に対して抵抗する力がある.この高分子の絡まりあいは急激な変形に対しては弾性挙動を示すが,応力がかかり続けると絡
まりあいがほどけ,時間と共に応力が緩和する.これが応力緩和の発現機構である.
このように,粘着剤の粘弾性は高分子材料であることに起因する.
また, ‑事由引っ張り測定においても粘着剤の粘弾性を見ることができる. Fig.2.10に同じ粘着剤で引っ張り速度を変えた時のSS曲線を示す.図中の曲線の右に書かれている数字は引っ張り速度である.引っ張り速度が速い場合は高分子鎖の絡まりあいがほどけない間に引っ張られるため応力は高くなるが,引っ張り速度が遅い場合は引っ張っている間に応力が緩和してしまい値は低くなる.このように,同じ粘着剤であっても引張り速度が変わるとSS曲線は異なる.つまり,粘着剤の‑軸伸張変
形には時間依存性のある粘弾性があることが判る.

0.7

貫o・6
i 0.5
む
a

・a
E=

oa

≡

萱[*:
0.2

100%

200%

300%

400%

500%

600%

700%

800%

900%

1000%

nominal strain

Fig. 2.10: SS behavior at verious tensile speed

このように,粘着剤の変形は複雑であり,単純な弾性体としては取り扱うことが
できない.
14

2.2.3 ピール試験での変形
粘着テープは多種多様な場面で用いられていることは上述した通りである.そのため,要求される必要特性も多様である.シミュレーションによる解析においても様々な実用評価を解析する必要があるが,その代表的な例としてピール試験による引張強さを示す.

Fig. 2.ll: Peel measurement equipment

ピール試験について簡単に説明する.粘着テープを被着体に貼りあわせ,その端部をFig.2.11に示すような装置を用い,引き剥がす時に発生する抵抗力を測定するものである.用途によって,引っ張る速度,角度,雰囲気温度を変えて測定する.図
は粘着テープと被着体とのなす角度が900であることから, 900ピール試験と呼ば
オI'J.

粘着テープがピール引張り試験により被着体から引き剥がされている実際の様子
をFig.2.12に示す.写真からも,外力が加わると粘着テープの粘着斉帽βが大きく変形していることが判る.
このように,粘着テープを引き剥がすために与えられた外力は,粘着剤が変形することにより分散し,被着体と粘着剤のより広い界面において剥離しようとする外力に抵抗することとなる.つまり,粘着斉摘号大きく変形することでより強い力を加えなければ剥がれないこととなる.このことからも判るように,粘着剤の変形する様子を詳しく解析することは粘着剤の設計に対して極めて重要なことである.

15

Fig. 2.12: Con丘guration in peel measurement

2.3 粘着剤の解析
これまでに,粘着剤を数値解析することの重要性を示した.それと同時に,粘着剤は複雑な変形挙動を示すこと,及びその変形は大変形を生じることを示した.このようなバルクの変形挙動を解析するには有限要素法による解析は有効である.しかし,粘着剤の複雑な変形挙動を記述できる適切な材料モデルが存在しない.また,大変形するため固体の有限要素法解析で通常用いられるLagrange型解法では計算格子が破綻してしまい解析することが難しい.

Fig. 2.13: Lagrangian computational result in peel

Fig.2.13にLagrange型解法でピール引っ張り試験を解析した例を示す.図中,下部の変形していない部分は被着体であり本解析ではSUSの鋼板を用いている.また,16

上部の変形していない部分は基材であり本解析ではPETを用いている.中央部の変形しているところが粘着剤であり,大きく変形していることが判る.また要素に着目すると',最初は長方形に分割されていた要素が解析の進行にともない形状が変形し,図では極めて歪な形状になっていることが判る.この時点で格子の形状破綻から計算は終了した.このように, Lagrange型解法を用いた解析では粘着剤の変形を解析することは難しい.
次に,粘着剤の解析に関する今までの検討について概観してみる.

2.3.1 粘着特性の解析
まず,剥離強さの研究を概観する.剥離強さの先駆的解析はRivl叫20],畑21],
およびDeryaginとKrotova[22]によりなされた. Rivlinは1944年に剥離過程にエネ
ルギー保存則を適用し,剥離強さが接着仕事と剥離角度の関数(式(2.10))であり,単位面積当りの仕事の次元を有することを示した.
Fp‑wa
1 ‑cosu(2.10)

ここで, Fpは剥離強さがなすエネルギー, waは接着仕事,のま剥離角度である.
畑は1947年に,曲げモーメントの釣り合い及び剥離仕事の平衡の計算方法から剥離強さがなすエネルギーを定式化している.曲げモーメント法からは接着仕事と剥離角度のみならず,表面基材の曲率半径と厚さの関数として剥離強さがなすエネル
ギーを表した(式(2.ll)).
wa

Fp‑(2.ll)
8p+cz ‑ cos6

ここで, ββは剥離時の曲率半径であり, czは表面基材の厚さである.剥離時の曲率半径が厚みに比べて十分大きな場合や,厚みを無視した取り扱いをすることにより, Rivlinの式(2.10)が得られる.
また,剥離強さの式は弾性力学を用いた解析からも得られ, Bikerma可23], Gardon囲,
Kealble[25]らによって定式化されている囲.三木は連続体力学から剥離理論の定式化を試み,剥離に関してレオロジー理論と界面化学理論を結びつける関係式を得ている[27].さらに, Gentは破壊の観点から粘着強さを論じている囲.

17

実際の実験結果と粘弾性試験から得られた物性との相関に関して, Daleらは室温の特性であっても高温環境下で得られた物性の方が室温で得られた物性よりもより実験と相関があったことを示している[29].
次に,粘着剤の粘弾性に深く関係があるタックの研究に関して概観する.タック
は粘着剤のもつ"ベトベト感"を数値的に評価したものであり,プローブタック試
敬,ループタック試験,ボールタック試験など試験方法は種々存在するが,詳細は参考文献に譲る[30].
タック発現の先駆的な研究はDahlquistによって行われ,タック発現のためのDahlquist
Criteriaを提唱した囲. ChangはDahlquist Criteriaを発展させ,動的粘弾性試験から得られる貯蔵弾性率と損失弾性率を用いて異なる粘着剤におけるViscoelastic
Windowと呼ばれるタック発現可能な領域を提唱した[31]. Davidは,貯蔵弾性率は粘着剤が被着体にくっ付く時に,損失弾性率は粘着剤が被着体から剥がされる時に影響することを論じている[32].

2.3.2 材料モデル
本研究ではマクロスケールを取り扱うため,まずマクロスケールでの材料モデル
について概観する.
粘着剤に対して粘弾性モデルを適応させた研究例は多い.それらでは, 2つもしくはそれ以上の複数のMaxwell要素を並列に並べたモデルを用いている.畑[33],水町囲[35]はこのモデルを用いて,粘着剤に一定応力が作用する場合の時間と共に粘着剤のひずみが増加するクリープ現象を論じている.
粘着剤は粘弾性以外にゴム弾性を有することは上述した.ゴム弾性を記述する材料モデルとしては超弾性モデルがある.超弾性モデルにはneo‑Hookeanモデルに始まり, Mooney‑Rivlinモデル囲, Yeohモデルなどがある.大変形下におけるゴムの挙動を精度良く再現できるモデルとして,近年ではOgdenモデル岡がよく用いられている.これら超弾性モデルを用いてSimo囲, Sussman囲,山田囲らはゴムなどの非圧縮性材料の解析を行っている.
次に,分子の動きというミクロスケールから成り立つ材料モデルについて概観す
る.多くの質量点をバネで結合したモデルで分子を表すRouseモデルが提唱され,このRouseモデルに質量点での抵抗を設け分子の通り抜けを排除したZimmモデル
18

がある.希薄溶液中での粘弾性挙動はZimmモデルに良く一致する. de Gennsが提唱し土井, Edwardsらが発展させたレプデーションモデル[叫では,粘着剤のような高分子からなる高粘度状態での高分子の挙動をよく記述できる.
このような分子の挙動から考えられたモデルにArrudaらが提唱した材料モデルがある[42].富田らはそのモデルをゴムの変形解析に適応させた研究を行っている[43][44].

2.3.3 シミュレーションに関する既存の研究
粘着剤の設計をする際には,粘着剤中の応力がどのように分布しているのか.また,様々に加えられる外力や変形に対してどのような応力分布を生じ,粘着剤がどのように変形するかについて知ることは非常に重要である.このような要求に対しては,有限要素法などによる数値解析をすることが望ましい.
粘着剤の変形を有限要素法で解析するには,粘着剤の応力‑ひずみ関係が複雑なこと,及び粘着剤が容易に大変形することによる難しさがある.
三木らは構成方程式としてLodgeモデル国を出発点とし,粘着剤の特性の1つである粘弾性モデルを導入して, 、大変形に適応可能な構成方程式を示している囲.
Duncanは粘弾性の物性に関してProny級数を用いて定め,その値を用いて粘弾性体のクリープ変形解析について有限要素法による解析を行っている[47].
朱峰は粘着剤の変形を3つの領域に分け,その各々に異なる解析手法を適応する
ことで,粘着剤の剥離を有限要素法により解析している囲.この手法は粘着剤のピール剥離試験の解析に特化したものとなっているが, Kealbleが実験で求めた粘着
テープの剥離時における剥離界面の応力分布国を定性的に表現できていることは興味深い.
山口らは粘着剤を引き剥がす際に生じるキャビテーションについて数値解析の観点から研究しており,プローブと被着体に挟まれた粘着剤を引き剥がす過程において,応力曲線やキャビティ生成の場所に関して実験と定性的に一致している結果を得ている囲.
近年, OCTA剛などの出現によって,分子動力学を用いたメソスケールのシミュ
レーションが進展しており,この技術を粘着剤に適応している例も増えている[叫.

19

第3章

Euler型有限要素法

本研究では,粘着剤が変形する際に生じる大変形を解析するために,固体解析に通常用いられるLagrange型有限要素法ではなく,大変形しても計算格子が歪まない
Euler型有限要素法を用いる.
本章では, Euler型有限要素法について説明する.まず,有限要素法の基礎方程式
である保存則について説明する.次に, Euler記述の方程式を解くための手法である
Operator Split法について説明する. Operator Split法は,非移流ステップ(Lagrangeステップ)と移流ステップ(Eulerステップ)の2つのステップに分割して解く方法である.最初に,非移流ステップについて説明する.本研究では,非移流ステップには陽解法を採用している.その後,移流ステップについて説明する.本章では,移流方程式の解法として, 2次精度を有するMUSCL法囲囲図を説明する.最後に,界面補足法である密度関数法[55][56]について説明し,時間増分制御につV)て述べることにする.

3.1 保存則
物理学は保存則を基礎として構成されており,連続体力学においても保存則は重要な役割を果たす.連続体力学における保存則には,質量保存の法則,運動量保存の法則,角運動量保存の法則,エネルギー保存の法則がある.本節では,特に固体解析において基礎方程式として用いる,質量保存の法則と運動方程式に対応する運動量保存の法則について説明する[57H64.

3.1.1 質量保存の法則
物質の質量mは,質量密度をp,物質の占める領域をVとして,
m‑ /pdv

20(3.1)

により与えられる.質量保存の法則は,質量mが時間に依存せず,変形後も一定,
すなわち
(3.2)

が成立することを述べている. D/DHま,物質時間導関数である.変形前の体積をdV,変形後の体積をdv,変形勾配テンソルをFとする.このとき,これらの間には,
dv ‑ (detF)dV(3‑3)

なる関係が成り立つ.ここで, detFは体積変化率を表し,
detF‑J(3.4)とする.式(3.3),式(3.4)より,式(3.1)は

m‑pJdV(3.5)と書き換えられる.これを式(3.2)に代入すると,
f(慧J+pDJ
Dt,dV‑0/(完J+JptiL)dV‑0/(霊ptiLJdV‑0/(完+p∇蝣v)dv‑O(3.6)となる.ただし, Vは物質点の速度ベクトルを表す.上式は物質の任意の一部分についても成り立つことから,

完+p∇・v‑ 0(3.7)

を得る.式(3.7)は連続の式と呼ばれる.非圧縮性物質の場合,質量密度βの物質時間導関数は零となることより,連続の式は次式になる.
∇v‑0

21

(3.8)

3.1.2 運動量保存の法則
物質に作用する力には,物体力bと表面力tがある.ただし, bは単位質量当たりの物体力, tは単位面積当たりの表面力とする.運動量保存の法則により,物質全体
における物体力と表面力の和と運動量の物質時間導関数は次式のように等値される.

芸 pvdv)‑ /pbdv+/ids(3.9)

式(3.9)において/ pvdiは物質全体の運動量の物質時間導関数Ivpbdtは物質全体における物体力の総和Jstdsは物質全体における表面力の総和を意味している・
式(3.9)はEulerの第1運動法則と呼ばれる.ここで,式(3.9)の左辺は以下のように変形される.

芸 pvdvI‑芸 pvJdVdV‑/(慧vJ+p芸J+pv霊)
pfuv+ ( />7r J+pv{∇ v)JβP
p芸dv + ( j元ev+.!‑ p芸dv + V (芸+p∇・vldv‑ /p芸dv
ββ∫.

∫.

(3.10)

pv(∇・V).

また,式(3.9)の右辺第二項は以下のように変形される.
tds‑‑T蝣nds(3.ll)

ここで, UはCauchy応力, nは任意面S上の外向単位法線ベクトルであり, Cauchyの公式とGaussの発散定理を用いた. Cauchyの公式とは,任意面の表面力(応力ベクトル) tが,任意面の外向単位法線ベクトルnのCauchy応力テンソルJによる#.'こ;射i'kt=an(3.12)

により求められることを述べている.なお,角運動量保存の法則により,
‑

22

a

(3.13)

なる関係式が成り立つことが要請される.つまりCauchy応力は対称テンソルである.
以上の式(3.10)(3.ll)より, Eulerの第1運動法則(3.9)は以下のように書き換えられる.

/p芸dv‑ (pb+∇・er)dv/p(芸‑b十<r)dv‑0(3.14)(3.15)

上式は物質の任意の一部分についても成り立つことから,

p芸‑∇ <r+pb
pa‑∇

pb

を得る.ここで, αは物質点の加速度ベクトルである.式(3.17)は, Cauchyの第1運動法則または平衡方程式と呼ばれる.

3.2

0perator Split法

弾塑性材料などの大変形固体解析で用いられる構成方程式の応力評価は,各物質点に追従した点で定義されるひずみ速度が用いられる.しかしながら, Euler表示での平衡方程式では速度の空間時間導関数が用いられている.よってEuler表示の平衡方程式を解くことによって求められたひずみ速度を,そのまま構成方程式で用いることはできない.このような問題点を解決するために,ここではEuler表示の平衡方程式に対してOperator split法[65]を用いる.
Euler記述の平衡方程式は次式で与えられる.

)可‑∇ cr+pbp(雷+{v‑∇)V)(3.18)

Fig.3.1は時刻nからn+ 1へと変形を進める際のOperator split法の概念図である. Operator split法では,式(3.18)を以下の2つの部分に分割する.

p(雷)ム/st.'B

‑∇

+pv

23

cr+pb

∇

<r‑0

(3.19)(3.20)

■
(.
し
+
蝣

帆

鵜
A
ソ

^ I V

s

午
!蝣

磨 妻
S+
・
I
JYガ
;=Yl
, 千
、
ニ
H
r

Non‑advective

蝣

謬
+蔓

Advective step

JI

Jl+I

Fig. 3.1: Concept of Operator split method

式(3.19)は外力項を含んだLagrangeステッフO,式(3.20)は移流項を含んだEuler
ステップである.また,これらの時間変化率の関係は次式で表される.

霊‑(訂+(霊E(3.21)

式(3.19)の上添え字Lで示された速度め時間変化率は,各物質点に追従した点で観測されるLagrange的な変形速度の時間変化率を示している.このLagrangeステッ
プを解くことによって求められたひずみ速度を用いて,構成方程式より応力を算出する.一方,式(3.20)で示された速度の時間変化率は, Euler的な変形速度を,壁間に固定されたメッシュへ投影する際の変化率を示している.

3.3 非移流ステップ
3.3.1 平衡方程式の離散化
Lagrange記述の平衡方程式(3.17)は,
paL‑∇

24

cr+pb

(3.22)

である.次に,平衡方程式から仮想仕事式を導く.平衡方程式(3.22)の両辺に任意の仮想変位6uを内積する.
pa,L蝣5u‑∇(3.23)

Su+pb蝣Su

ここで,
∇

5u)‑

∇

a)蝣5u+cr:(Su㊧∇)(3.24)なる関係が成り立っことを利用すると,式(3.23)は以下のように変形される.
pai,蝣Su‑∇・<7‑∂u)‑cr:(5u㊥∇pa,L蝣Su+a:(Su㊧∇)‑∇pb蝣Su(<x・Su)+pb蝣Su

上式を領域0で積分すると,次式を得る.
/paL蝣8udV+/<x:(8u'nJn㊤∇)dv‑I
Jn∇Su)dV+/pb‑SudV(3.27)
ここで,右辺第1項に対してGaussの発散定理を適用し以下の式を得る.
Jn
paL蝣SudV+/<r:(8u∇)dV‑ft‑sudS+fpb‑SudV(3.28)
JvtJn
ここでSは物体の表面領域であり,自ま表面カベクトルである.式(3.28)は仮想仕事式と呼ばれる.
仮想仕事式(3.28)を,有限要素法により離散化する.後述するアイソパラメトリッ
ク要素を用いて離散化を行うと,結果的に次の離散化方程式が導出される.

Ma+Fint‑Fe蝣ext(3.29)

ここで, Mは対角化された集中質量行列, aは加速度ベクトルであり intと 蝣ext
は,それぞれ内力および外力ベクトルである.

3.3.2 動的陽解法
本研究ではLagrangeステップにおいて動的陽解法囲[67]を採用する.動的陽解法において,物質点の位置ベクトルを移動させることなく,応力などの経路依存のある変数の更新が可能である.よって通常の固体解析のためのLagrange記述における動的陽解法を,そのまま用いもことが可能である.
25

今,時刻tnでの力学挙動が既知であるとして, △ま秒後の時刻tn+lでの力学挙動を求めることを考える.ここで,一般的な関数fに対してテイラー展開を行う.ただし,時刻tnと時刻tn+1の中間時刻を時刻tn+芸とする.
nfn+
r‑fL悪(芸△恒7去
1∂2fn+/'/‑'(芸△tf ‑

(3.30)(芸△t?+‑(3.31)

lf去
1∂2fn+

fn+l‑fn+¥+悪(芸△t)‑2∂ま2△t2のオーダーの項まで考慮して,式(3.30)から式(3.31)を引くと以下のように
K 'J.

r+l ‑ fn +悪(△t)(3.32)

式(3.32)は中央差分法といわれる手法[52]であり,これを適用すると時刻un+1における位置ベクトルx71+1は,時刻tn+吉での速度ベクトルvn+吉により,次のように評
価できる.
xn+l = xn +vn+^△(3.33)

また,時刻まn+吉での速度ベクトルvn+圭は,時刻tnでの加速度ベクトルanを用いて次のようになる.
vn+去‑ vn一言+an△(3.34)

ここで,時刻tnでの加速度ベクトルanは,式(3.29)より次のように求められる.

a'‑MJv‑1ext iint)(3.35)

以上より,動的陽解法の流れを以下に示す.
1.時刻tnでの加速度ベクトルを求める.
2.時刻tn+壬での速度ベクトルを求める.
3.時刻tn+lでの位置ベクトルを求める.
4.時刻tnから時刻まn+1までの応力積分を行い,時刻tn+lでの応力を求める.
この1から4の流れを繰り返して解析を実行する.
陽解法では反復計算が不要であるため, Newton‑Raphson法囲等による反復計
算を用いた陰解法(静的問題の弧長法や動的問題のNewmark‑β法などでありがちな,26

計算途中で収束が得られずに解析が停止するブレークダウンの頻度が少ない.また,
計算スキームが非常に簡便である.例えば,陰解法において健全な収束を得るための接線剛性マトリックスの導出に際しては,非常に煩雑な理論展開が必要となるが,陽解法ではその接線剛性マトリックスが不要である.そのため,解析で必要とされるメモリが陰解法に比べて少ない.さらには,解析の並列処理を考えた場合,内力および外力ベクトルのみを分割すればよいので,並列処理の実行が容易などの利点がある.しかしながら,陽解法で安定した解を得るためには,時間ステップサイズ
をある値以下にしなくてはならない条件, Courant条件が課せられる囲囲.そのために時間スケールの長い問題は膨大な時間ステップが必要となる.

3.3.3

アイソパラメトリック要素

本節では,本研究で用いているアイソパラメトリック要素[66][67]について説明
する.
問題を記述する物理座標系x,y,z)とは別に自然座標系&v,0を導入し,図312のように(LvX)をパラメータとする立方体領域から変換されて,任意の8節点六面体要素が得られる場合を考える.この時,自然座標から実際の物理座標系におけるx, y, zへの変換を表す式は次のようになる.

(3.36)
・‑{^(」,77,0^(^,0・・・N8(Z,V,()y‑{^(」,77,0 ^(e両) ‑ Ns(」,相〉

27(3.37)

・‑{wi(」,相N2(tv,C) Ns(Z,相〉
ここで(3.38)

Na(Z,V,0は形状関数と呼ばれ,次のように表される.
・1‑芸(1‑0(1‑サ7)(1‑C)

・2‑芸(n‑e)(i‑^)(i‑C)
・3‑芸(1+0(1十サ?)(!‑0

・4‑芸a‑oa+^a‑o(3.39)

・5‑芸(1‑」)(1‑77)(1+C)
・6‑芸(i+0(ト*7)(1+C)
・7‑芸a+oa+^a+o
・8‑芸(l‑0(l+サ7)(l+0
式(3.39)において&v,0は要素領域を{ォ蝣刷

‑1≦E≦1, ‑1≦り≦1, ‑1≦C≦(3.40)とする.また,式(3.36)(3.37)(3.38) (3.39)において,右下添え字の1‑8は, Fig.3.2に示すアイソパラメトリック要素の節点番号である.
本研究で用いるアイソパラメトリック要素は,これらの座標変換で用いた形状関数を物理量の内挿関数として用いる要素である.したがって,任意の形状をした8節点六面体要素の節点変位を用いて,各変位成分を次のように表すことができる.
8

u ps Nu ‑ ∑Na&両)iiα7&¥(3.41)

前述のように,本研究では,これらの近似した変位を用いて仮想仕事の原理を離散化し,離散化方程式(3.29)を得る.

28

8

7

(
5

‑

一
一
一 {ト
午I
6
77
+ L ‑

1‑'
0

rE

4

3
‑ 1

1

2

Fig. 3.2: 3D 8‑node isoparametric element

3.4 移流ステップ
移流ステップでは,密度関数￠,速度V,応力Uに関する移流方程式を解く.これらについて以下で説明する.

3.4.1 密度関数
既述のように,質量保存の原理(3.2)は,
Dm
‑=O
Dt

(3.42)であった.後述する密度関数を￠とし,有限要素の体積をVele 物体の質量をpとすれば,ひとつの有限要素に対する質量保存の原理は,D (pVeeleSDt=0(3.43)となる.ここで, pとVeleは時間に対して一定であるため,

慧‑o(3.44)

雷+V・霊‑o(3.45)となる.本研究では,密度関数に関する移流方程式(3.45)を解くことにより物体の境界面を表現している.
29

3.4.2 速度
Lagrange記述の平衡方程式(3.17)は,次式のように表現することもできる.

p芸‑∇ cr+pb(3.46)

ここで, Vは速度である.さらに,これをEuler記述に書き換えると,∂V

諒+V霊‑去∇ cr+b(3.47)となる.既述のように, Operator Split法ではEuler記述の方程式(3.47)を以下のように分割する.

塑=旦∇ a+b(3.48)

∂t p
∂V

川

∂V

r蝣一石

0

(3.49)

式(3.48)を解くのが非移流ステップ(Lagragianステップ)であり,式(3.48)を有限要素法により離散化した方程式が前述の方程式(3.29)である.一方,式(3.49)は速度Vに関する移流方程式であり,これを移流ステップ(Eulerステップ)で計算する.

3.5
3.5.1

VOF法による界面捕捉
VOF法

voF法は, VOF関数￠を要素毎に定義し,その値によってその要素にどの程度の割合で物質が存在しているかを表す方法である.すなわち, ￠‑0であれば要素内にその物体は存在しない,または他の物体が存在することを意味し, ￠‑1であれば要素はその物体で満たされていることを示す.また, o<￠<1であればその要素に境界面が存在することになる.
voF法では,以下のVOF関数の移流方程式を導入し,これを解くことでVOFを移流させる.

慧‑霊+V・∇桓0(3.50)

ここで,∇・(V￠‑V

∇'+

30

∇・V

(3.51)

の関係を考慮すると,式(3.50)は以下のようになる.

霊+∇・(V￠ト￠∇蝣v‑0(3.52)

非圧縮性を仮定すると,非圧縮性の連続の式が成り立つ.したがって,上式の左辺第3項は0となり,式(3.50)は以下の式に置き換えることができる.

霊+∇・(V￠ ‑o(3.53)

上式の形に変えたのは,次節で説明する発散型スキームのMUSCL法を適用するためである.

3.5.2 2次精度風上差分法
2次精度風上差分法として, MUSCL法を用いる. VOFのなまりが拡大しすぎると,物体の変形状態を捉えることが困難になる.一般的には式(3.53)を解く際に,精度保持のためには1次精度では不十分であり,最低2次精度は必要とされると考えられている囲.したがって,式(3.53)の解法として2次精度のMUSCL(Monotone
Upwind Schemes for Conservation Laws)法囲図を用いる. MUSCL法では,要素境界における両側の物理量の数値流束を用いてその勾配を評価する.このときに要素両側の物理量として,両境界に接する要素の値をそのまま使うのではなく,さらに隣りの要素での値を用いて定めた,要素内部の物理量の分布関数によって求める.このとき,制限関数を用いることで,後で述べるTVD条件と呼ばれる条件を満たすように調整する.こうすることで数値振動の発生などの不具合を抑え,計算を安定に進めることができる.

TVD条件
2次精度以上の移流方程式の解法を用いる場合,計算が進むにつれてVOFの数値
振動が大きくなることがある.この数値振動は物理量の種類によって非現実的な計算結果をもたらす.このように,数値振動は数値計算にとって致命的になることがあるので,移流方程式の解法には単調性を維持できるスキームが必要となる.単調mであるという条件を緩めて,全変動(total variation)が増加しないという概念が提案されている囲.全変動は次のように定義される.

TV(の‑∑ Di+1‑q)A(3.54)

7,

この全変動とは自由・移動境界面の凹凸具合を示すものである.この全変動が増加しないとする条件をTVD(Total Variation Diminishing)条件榊と言い,以下のように表される.

TV{(pn+l) < TV(の(3.55)

つまり,この条件を満たす解法は単調性を維持できる囲. MUSCL法では,後で述べる制限関数を導入してTVD条件を満足させている.

移流方程式の離散化式(3.53)を時間差分近似し,以下のように表す.
kn+l ‑△t∇vォbn)(3.56)

上式を,成分を用いて以下のように表す.

un+1‑￠m‑△t(笠+笠十等)(3.57)

さらに x,y,zの各方向の勾配に分けて,以下のように計算する.

・*‑￠n‑△t(笠)・**‑￠*‑△t(箸)(3.58)

kn+l ‑ ￠** ‑△t (笠)すなわち, 3:方向の勾配を計算した4,*を用いてy方向の勾配を計算して</>**を求め,
これをZ方向の勾配の計算に用いる.
式(3.57)を要素境界におけるVOFの数値流束の勾配で差分近似すると,以下の32

ようになる.

畔.+1‑

た‑△t

面i.去,j,k4>i+去,3,k 面i寸i,k(f>i寸,.,良

△T

‑△i

否i,j+去,k<f>i,j・3+喜,k‑否i,3‑圭,た0ォ一基,Lc△t

訂i,i,k+S<Pi,3,良+去‑両軸た‑主転た一主△Z)(3.59)

ここで,面,有声は速度の平均値であり,要素面中心の速度である.また, 4,は要素面中心におけるVOF値であり,これは次節で説明するように制限関数を用いて, 2次
精度的に内挿される.これらの位置関係をFig.3.3に示しておく(2次元).

Fig. 3.3: Advection of density function

制限関数化
本研究では, MUSCL法における制限関数としてSlope Limiterと呼ばれる制限関数を用いる.要素境界の物理量を内挿する際に, 2次精度のスキームでは要素境界に不具合の原因となる極大が発生してしまうことがある. SlopeLimiterは,この極大の発生を抑えてTVD条件を満たす働きをする.
要素境界の物理量を内挿するために, Fig.3.4に示すように,各要素に傾き(Slope)を定義する.なお, MUSCL法では多次元空間のそれぞれの方向で1次元的移流計
算を行う.したがって,説明は1次元で行うことにする.

33

Fig. 3.4: Slope Limiter

定義したSlopeを用いて,要素境界の物理量が以下のように計算される.

転

si‑i筈

勘‑￠i‑SiAx(3.60)

娠‑恒轄
職 j+1‑s,.i+1害

ここで,上添え字のRとLはFig.3.4で右から内挿した場合と,左から内挿した場
合との違いを表しており,その点での速度の向きによって使い分け,風上化を施す.
Slopeは以下の式により与えられる.

S ‑ ‑ {sgn(s/) + sgn(sr)} min(￨s/￨, ¥sc冊
ここで. sgn(F)はFの符号を求める関数である.また(3.61)

sl,sc,srは要素iにおいては以下のように計算される.
sl=￠i ‑ <pi̲1

△z
i+1

,<(・

,SI

=

2△z
￠i+1
=
△3;

(3.62)

式(3.61)より得られたSを用いて,式(3.60)により要素境界の物理量を求める.風
上化を施すために,速度が正の場合は左から内挿したものを使い,負の場合は右から内挿したものを使い,数値流束子を計算する.このとき,速度の向きを判別しなくても良いように,以下のようにまとめる.
7, ...'二ト̲/,':I‑万...' ,.'・・,

fi‑喜‑

34

(3.63)

if. !言I一了∴'・蝣) I万両一上∴'"I'
fi.去‑(3.64)

この数値流束の勾配を評価することにより,移流計算を行う.この時,解析領域の端部要素のVOFについて計算するためにゴーストメッシュを設置している. VOFの移流は2次精度で計算するため, Fig.3.5に示すようにゴーストメッシュは2つあり,その両方に端部要素のVOFと同じ値を入れている.

Fig. 3.5: Ghost mesh for MUSCL method

i+1

i+2

Fig. 3.6: Slope Limiter

最後に,式(3.61)の働きについて説明する. VOFの分布がFig.3.6のような場合を考える・このような場合,本来は最も小さな傾きであるscを勾配として,鴎主を計算する・しかし,そうすると￠i>転1>呪となり・要素境界に物理的には不合理な極大ができてしまうことになる.これはTVD条件に反することであり,敬値振動の発生につながる.したがって,これを防ぐために式(3.61)では,このような状況になった場合は, s‑oとなるように計算されるようになっている. 5‑0と35

いうことは,一次風上差分に落としたということになり,鶴‑・>i+lとなる・そう
すると, ￠i>亀≧￠i+1となり,単調性が維持できる・

3.6 時間増分制御
本節では, Euler型解法で用いられる時間安定条件について説明する.
動的陽解法の時間積分において,時間増分△舌を,ある限界値を超えて大きくとると,解が桁違いに不合理な値をとり発散に至る場合がある.
このような不安定現象が△tを小さくとることで回避される場合,その時間積分法は条件付安定であると呼ぶ.
一方, △tの値のいかんによらず安定な解が得られる場合,その手法は無条件安定であると呼ぶ.
本研究で用いている中央差分法に基づく陽解法は,条件付安定である.この場合,△tは要素内を応力波が伝播する時間より小さくならなければならなし子.すなわち,△t <車型(3.65)rォ

なる条件(Courant条件と呼ばれる)が課せられる.ここでLmmは解析対象を構成する要素の中で最小の接点間距離である.また, Cは応力波が伝播する速度であり,式(3.66)で与えられる.

C‑据(3.66)

ただし,別まヤング率, βは密度である.
さらに,移流方程式を解く際の時間増分を△t/とすると, △ま′は,要素内を物体が速度Vで伝播する時間より小さくならなければならない[52].すなわち,△t′<車型回(3.67)

なるCourant条件が課せられる.
本研究では,式(3.65)の条件と式(3.67)の条件を満足するために,
△t<a,

Lr,

max(c,柄)(3.68)

なるCourant条件を課している.ここで, CuはCourant係数と呼ばれo < Cu < tt/4を満たす値を用いてスケーリングを行う.
36

第4章 粘着剤の構成方程式
本章では,粘着剤の変形挙動を記述する構成方程式について説明する.粘着剤の変形を記述するには,粘弾性とゴム弾性の両方を記述する必要がある.まず,粘弾性を記述する構成方程式を構築する.その後,得られた構成方程式についてゴム弾性を記述できるよう拡張することについて示す.

4.1 微小変形理論における等方性弾性体の構成方程式
微小変形理論における等方性弾性体の構成方程式は, Cauchy応力テンソルU,微小ひずみテンソルE,弾性係数テンソルCを用いて,(4.1)⇔<T‑2/i」入tie)I(4.2)と表現できる岡[叫囲.ただし,この構成方程式i.1)は,基準枠(観測者)が変
われば,この方程式も変わるため,物質客観性の原理を満たしていない.物質客観性の原理については, 2.3.3節で詳説している.式(4.1)において, FLと人は等方怪
弾性体のラーメの定数と呼ばれ,等方性弾性体の各種定数には,表4.1に示す関係がある[70].
また,微小ひずみテンソルEを用いて,偏差ひずみテンソルは

E′

=e‑

‑(tre)/

37(4.3)

Table4.1: Relations between various constants of isotropic elastic material

定数の呼称

入,pによる表示

入 ラーメの定数

入Kv,Gによる表示/(.一一・('E,のこよる表示(1+1/ 1‑2!/

H‑G
せん断弾性係数

E

〃

丁丁丁

Kv

Kv体積弾性係数 入十号〃
E

ヤング率

9KvG
3Kv+G

M3人+2a*S9Z

E
3(1‑2;/)

E

3Kv ‑2G
2(3∬+C)〃 ポアソン比で定義される.この式(4.3)を用いれば,式(4.2)は,cr‑2G e'+l(tre)I+A tre /‑2Gs'+呈G+X)(tre)J
=

2Ge′+

(4.4)

tve)I

(呈G+Kv

・,‑ 2Ge'+ Kv(tre)Iとなる.この式(4.4)の右辺第一項は偏差応力テンソル,第二項は圧力項を意味する
から,
J′

≡

2Ge'

p ≡ ‑Kv{skk)と定義すれば,式(4.4)は
cr

‑

(T′

‑pi

(4.7)

と変形できる.すなわち,微小変形理論における等方性弾性体の構成方程式(4.1)は,何の仮定をおくことなく,式(4.7)に変形できることから,式(4.1)と式(4.7)は等価て'.¥t"ら38なお,以上の構成方程式において,ラーメの定数Åまたは体積弾性係数Kvは,秦4.1より,ポアソン比L/とヤング率Eを用いて表示すれば,
Eu
Å=

(4.8)
(1+〟 1‑2〃)

Kv‑

(4.9)
3(1‑2〃)となる.よって,非圧縮性物質すなわちポアソン比が1/2になる場合は,これらの定数は∞になり,以上の構成方程式は適用できない.
粘着剤は非圧縮に近い性質を示すが,完全な非圧縮を示すものではない.よって粘着剤の特性を微圧縮性として構成方程式を考える.
よって,以下では微小変形理論における等方性弾性体の構成方程式(4.7)を用いて,微圧縮性粘弾性体の構成方程式を導くことにする.

4.2

一般化Maxwellモデル

粘着剤が有する材料物性の1つに粘弾性がある. 2.2.2節で示したように,粘着剤は顕著な粘弾性挙動を示す.粘弾性を再現する材料モデルとして一般化Maxwellモ
デルがある[叫[71]囲. Fig.4.1はその一般化Maxwellモデルの概念図である.一般化MaxwellモデルはHooke弾性体を示すスプリングとNewton流体を示すダッシュ
ポットを直列につなぎ,さらにそれを並列につなげたモデルである.

Fig. 4.1: Generalized Maxwell model

ここでEiは弾性係数, 7再ま粘性係数であり,緩和時間Tiは次式のように定義される.
(4.10)

なお上添字iはモデルの要素番号であり Fig.4.1ではn要素のモデルを意味している.

4.2.1

一般化Maxwellモデルの定式化

ここで1次元状態の一般化Maxwellモデルを, 3次元問題へ適用することを考える. Cauchy応力Uは偏差応力qlと圧力pにより以下のように分割できることを前節で述べた.
q‑J′‑pi(4.ll)

以下ではまず偏差応力,その後に圧力について述べる.いまスプリングを線形弾性体であるとして, i番目要素のスプリングのせん断弾性係数をG% 番目要素のス
プリングのひずみをE皇pとすると,スプリングの偏差応力J/iの構成方程式は次式のように示される.
J′

‑

2Gz

(4.12)

なおスプリング部分は非圧縮ではないがほぼ非圧縮と見なせる微圧縮性を仮定すろため,式(4.12)において本来はスプリングのひずみの偏差成分が用いられるのに対して通常のひずみを用いている.ここで,せん断弾性係数Giと弾性係数Eiは,ポアソン比L/iを用いて以下のような関係にある.
E%Gl=

2(1 〟)(4.13)

後述の非線形解析において扱いやすい速度形の構成方程式岡を導くために,式(4.12)の物質微分を考える.
l

>er′‑

Dt

‑ ^D¥v(4.14)

ここで, D芸pはi番目要素のスプリングのひずみ速度である・
式(4.14)の構成方程式は微小変形問題でのみ成立するものである.この構成方程式の大変形問題への拡張を考える.微小変形問題の構成方程式(4.14)において,応力速度晋は客観応力速度に置き換えられるものとする[66].物質客観性に関しては後述の4.2.2章に詳しく示す.客観応力速度については,その取り扱いが容易であるJau‑ann速度(筈{J)を用いることにする・その場合の大変形状態での構成方
程式は以下のようになる.

(筈(J)‑2&Dl,
40(4.15)

構成方程式(4.15)において,右辺のひずみ速度Dは物質点の座標の関数であり,物質点と共に移動する観測者から見たひずみ速度の定義に一致する. Cauchy応力のJaumann速度

Dcrn

( Dt

) (J)

に関しても移動する観測者から見た応力速度という意味合いを持つ.つまり大変形問題においては,応力速度とひずみ速度を両方とも客観性のあるものに置き換え,どれだけ物質がゆがんだとしても微小変形問題での構成関係が成立すると考える.よって式(4.15)内の応力速度とひずみ速度を関係つけるGiは微小変形問題におけるものと同一である.なおCauchy応力のJaumann速度(誓(J)と通常の応力速度(晋)との関係は次のようである・
rD(TH¥‑fDan‑W‑tT+ar‑W(4.16)(J)

ここでⅥ′はスピンであり,上記のひずみ速度DとスピンWは次のようにそれぞれ速度勾配の対称および反対称成分で定義される.

D‑芸〈(,v)+ ∇v(4.17)

W ‑芸〈(∇Vト(∇v¥T(4.18)

一方で,ダッシュポット部分は非圧縮性ニュートン流体であると仮定する.そしてそのせん断粘性係数をりiさらにひずみ速度をDbpとすると, i番目要素のダッシュ
ポットにおける偏差応力J′iの構成方程式は次式のようになる.
crn ‑ 2rfD昌(4.19)

ここで一般化MaxwellモデルのひずみEはi番目要素のひずみEiと同一であり,かつスプリング部分E皇pとダッシュポット部分'dpの和で表現できると仮定する.
e‑e甘

sp

'dp

(4.20)

さらにひずみ速度も同様の仮定が成り立つとする.
D ‑ Dl ‑D呈p+D昌(4.21)

そして式(4.15)(4.19)(4.21)を整理すると,一般化Maxwellモデルのi番目要素の速度形構成方程式は次のようになる.
(Dall'
vDt,(J)= 2GID‑tT'

2(1 vi]Ti

41(4.22)

偏差応力の速度形構成方程式が式(4.22)のように一般化Maxwellモデルの要素ごとに定義されるのに対して,圧力に関してはすべての要素を含めたものを定義する.
いま圧力については粘弾性の特性を考慮せずに線形弾性体を仮定する.その場合の圧力pは次式のように定義される.
p ‑ ‑Kvtre(4.23)

ここでKvは体積弾性係数であり,すべての要素の総和を意味する∑を用いて次のようになる.

Kv‑云
i‑l

El
(4.24)
3(1 ‑ 2i/*)

ここでnは一般化Maxwellモデルの要素数である.
式(4.23)の両辺の物質時間微分をとれば,‑ ‑KvtiD(4.25)となる.なお圧力pはスカラーでありスピンⅥ′の影響を受けないので,式(4.25)は大変形状態においても客観性の有無に関係なく同一である.
以上をまとめると Fig.4.1に示す一般化Maxwellモデルに対して,大変形状態にも対応した速度形構成方程式を次式のように得る.
J/〜

2(1 + i/)r'‑・丁‑ T(4.26)

4.2.2 物質客観性の原理
本節では,構成方程式の検証における議論を明確にするため,物質客観性の原理[57]囲囲岡について説明する.そのために,まず観測者の概念や客観性のある
スカラー,ベクトル, 2階のテンソルの定義を明確にする.その後,それらの議論に基づいて物質客観性の原理を説明することにする.

観測者の変換
観測者とは,空間内の事象の位置ベクトル㌶と時間tの組軽舟を観測する存在である.いま, 2人の観測者0と0*が,空間内の事象をそれぞれ軽,t}, {x*,t*}42

と観測したとする. 2人の観測者が運動する相対速度が光速に比べて十分に小さい場合はt‑t*と見なせるから,空間内の事象は観測者O: {x,t}
観測者O* : {㌶',*}と観測される.このとき,観測者0が観測する位置ベクトル3:iは,x*‑c(t)+Q(苦)・諾

(4.29)

なる式で,観測者0*の観測する位置ベクトルa:*へ変換することができる[57][74][59].
ここで. Q(t)は観測者0から観測者0*への回転を表す直交テンソル, C(申ま観測者
0*から観測した観測者0の位置ベクトルである.式(4.29)を観測者の変換と呼ぶ.

客観性のあるテンソルと観測不変テンソルの定義(客観性のあるスカラーの定義)
観測者0が観測する任意のスカラー(0階のテンソル)をα,観測者0*が観測す
る任意のスカラーをα*とする.これらについて
α*‑α(4.30)なる関係式が成り立つとき,すなわち,観測者によらずスカラーαの値が一定であるとき,そのスカラーは客観性があると呼ばれる.直感的に理解しやすい例をあげ
れば,物体の質量は観測者によらず一定の値になるから,質量は客観性のあるスカラーであるといえる.

(客観性のあるベクトルの定義)
また,観測者0が観測する任意のベクトル(1階のテンソル)をb,観測者0*が観測する任意のベクトルをb*とする.これらが,b*‑Q{t)蝣b

43(4.31)

なる式で関係づけられるとき,そのベクトルには客観性があると呼ばれる.式(4.31)はベクトルbの方向が観測者の変換によって不変であることを保証している.これについても,直感的に理解しやすい例を考えると,空間の2点間の結ぶベクトルは客観性のあるベクトルであるといえる.

(客観性のある2階のテンソルの定義)
次に,客観性のあるベクトルの観測者変換の式(4.31)から,客観性のある2階のテンソルの定義式を導くことにする.
いま,客観性のある任意のベクトルdとbを考える.ベクトルbからベクトルdへの線形変換,すなわち2階のテンソルをTとすれば,
d=T‑b(4.32)となる.
任意のベクトルdとbには客観性があるものとして考えているから,客観性のあるベクトルの観測者変換の式(4.31)から,

d*‑Qit)蝣db*‑Q{t)蝣bが成り立つ.これらの式(4.33) (4.34)を(4.32)に代入すると,
d*‑Q{t)‑T‑Q'it)‑b*となる.ただし, Q(t)は直交テンツルであるから, QJit)‑Q(4.35)

ま)となる[64]ことを利用した.よって,観測者0が観測する2階のテンソルTと観測者0*が観測する2階のテンソルT*の関係はT*‑Q(t) T Ql(t)(4.36)となる. Tは客観性のあるベクトルの線形変換であるから, Tにも客観性があると考えられる.よって, 2人の観測者から観測される2階のテンソルが式(4.36)で関係づけられるとき,その2階のテンソルは客観性のある2階のテンソルであると呼ばれる.
44

以上の議論からわかるように,客観性のあるベクトルの観測者変換の式(4.31)と客観性のある2階のテンソルの観測者変換の式(4.36)は,ベクトルbと2階のテン
ソルTに関連する方向が,観測者の変換によって不変であることを保証している[5].

(観測不変テンソルの定義)
観測者0が観測する2階のテンソルTと観測者0*が観測する2階のテンソルT*
について,

T*=T

(4.37)

が成り立つとき,このテンソルを観測不変テンソルという.すなわち,観測不変テ
ンソルとは,文字通り観測者の変換に対して不変なテンソルである.

客観性のあるテンソル・客観性のないテンソル・観測不変テンソル一覧表客観性のあるテンソル・客観性のないテンソル・観測不変テンソルについて,主たるものを以下の表4.2にまとめる岡囲[74][75].

45Table4.2: Objectivity of tensors

客観性のあるテンソル(T

Cauchy応力テンソル

V

左ストレッチテンソル

β

左Cauchy‑Green変形テンソル

A

Almansiひずみテンソルn

変形速度(ストレッチング)テンソル

(慧(J)
(雷(。)
¥Dt)(C)(宕(G)cauchy応力テンソルのJaumann速度
cauchy応力テンソルのOldroyd速度
cauchy応力テンソルのCotter‑Rivlin速度
cauchy応力テンソルのGreen‑Naghdi速度

客観性のないテンソル
̲F

変形勾配テンソル>r

微小ひずみテンソル/'IT

Cauchy応力テンソルの物質時間導関数

R

物質点の剛体回転を表す直交テンソル

L

速度勾配テンソルlγ

スピンテンソルfl

剛体スピンテンソル

観測不変テンソルu

右ストレッチングテンソル

C

右Cauchy‑Green変形テンソル

E

Green‑Lagrangeひずみテンソルs

第2Piola‑Kirchihoff応力テンソル

呈第2Piola‑Kirchihoff応力テンソルの物質時間導関数/ >1:Green‑Lagrangeひずみテンソルの物質時間導関数46

物質客観性の原理
ある物質が外力作用を受けると,その物質は特有の力学的挙動を示す.このように,各種の物質がそれぞれ特有の力学的挙動を示すのは,それぞれの物質を構成す
る分子や原子の構造などの相違が原因となっており,観測者には無関係である.
したがって,各種の物質の力学的挙動を記述する構成方程式は,観測者によらず同一の方程式にならねばならない.このことは物質客観性の原理と呼ばれるもので,構成方程式が満足しなければならない最も基本的な必要条件である.
よって,構成方程式が物質客観性の原理に矛盾しないためには,(1)構成方程式の両辺に客観性のあるテンソルを用いる.
(2)構成方程式の両辺に観測不変テンソルを用いる.
のいずれかである必要がある.
つまり,構成方程式に用いられるテンソルに「客観性がないといけない」というわけではない.物質客観性の原理という言葉の「客観性」は, 「観測者によらず構成
方程式が同一の方程式にならなければならない」という意味である.よって,構成
方程式の両辺に観測不変テンソルを用いた場合でも,物質客観性の原理を満たすことができる.一方,テンソルの客観性という言葉の「客観性」は式(4.36)を満足するという意味であって,物質客観性の原理の「客観性」とは意味が異なる.

47

4.3

般化Maxwellモデルの修正

前節では粘着剤の粘弾性挙動を記述するために,一般化Maxwellモデルを用いることを説明し,その構成方程式の定式化を行った.しかし,一般化Maxwellモデル
では粘着剤のもう1つの材料特性であるゴム弾性に起因する, ‑軸引っ張り試験において高ひずみ領域で応力が指数関数的に増加する挙動は再現できない.
これはFig.4.1の一般化Maxwellモデルにおけるスプリング部分の弾性係数を宝数としているためである.そこでこのスプリング部分の弾性係数を定数とする代わりに,モデルに作用するひずみEによって弾性係数は変化するものと修正する. Fig.4.2は一般化Maxwellモデルを修正した新たな材料モデルの概念図である.修正されたモデルではすべての要素の弾性係数EiがひずみEの関数となっているが,緩和時間
Tiは一定であるとする.よって緩和時間は修正前の一般化Maxwellモデルと同様に式(4.10)のように定義されるので,すべての要素の粘性係数77iもひずみEの関数となる.

Fig. 4.2: Modified generalized Maxwell model

4.3.1 弾性係数
ここでは,ひずみEとそれぞれの要素の弾性係数Eiとの関係について調べる.そのために一般的な粘着剤の応力緩和試験を行った.一方で一般化Maxwellモデルの応力緩和を示す理論式は式(4.38)のようになる.応力緩和の測定値と式(4.38)から得られる理論値を用いて最小二乗法による回帰分析を行い,すべての緩和時間に対応した弾性係数Eiを求めることができる.なおここではモデルの要素数は5とし48

て,その緩和時間は1.00, 1.00×10¥ 1.00×i(r, 1.00×10d, 1.00×106 の5種類でJ>'.).

6‑差」*exp一芸(4.38)

またそれぞれの緩和時間におけるひずみの変化による弾性係数の影響を調べるた
めに,初期ひずみを100, 200, 300, 400, 500, 600%の6種類設定して応力緩和の測定を行う.その結果をTable.4.3に示す.実際の測定結果は公称ひずみ一公称応力の関係で得られるが,表には真ひずみ‑真応力に変換した値を示している.
Table

4.3: Relaxation times and elastic moduli of each strain

結果的に,縦軸に求めた弾性係数そして横軸にひずみとする図を,緩和時間の種類分つまり5種類措くことができる.そのうちの1つとして, Fig.4.3に緩和時間が1.00×102 の場合の弾性係数とひずみの関係を示す.図中のプロットは100%から
600%までの6種類のひずみに対応する弾性係数である.

49

0.35
0.30
co
'oil
邑0.25
こ′l
戸
言0.20
CI
∈
0.15
胃
w0.10

0.0

0.5

,

1.5

2.0

2.5

strain [‑]Fig. 4.3: Relationship between elastic coefficients and strains in relaxation time
1.00×102

Fig.4.3においてプロット点は実験から求めた値であり,曲線は最小二乗法により指数関数でカーブフイテイングしたものである.図から判るように両者の整合性は非常に高いことが判る.また1.00×102s以外の他の4種類の緩和時間の場合においてもこの傾向は同様である.以上のことから本論文では,一般化Maxwellモデルの各要素の弾性係数Eiを,次式のようなひずみEの指数関数となるように修正する.
E%‑

Alexp(5号E)(4.39)

ここで, AiとBiは材料パラメータである.
このように弾性係数がひずみの増加と共に指数的に増加する理由を考えてみる.
物質の弾性には結晶弾性とゴム弾性があり,粘着剤の弾性はゴム弾性が支配的であることは2.2節で説明した通りである.ひずみを与えるとそれが小さい時は高分子鎖の絡まり合いが引き伸ばされ,エントロピーを増大させるべくゴム弾性が発現することを述べた.しかし,与えるひずみが大きくなると高分子鎖が伸びきった状態になり,さらにひずみを与え続けると,今まで殆ど生じることが無かった共有結合の距離や角度が変化し始める.つまり,結晶弾性の成分が生じることとなる.結晶
弾性はゴム弾性に比べて極めて弾性率が大きいため,少しのひずみの変化で大きな応力が発生することとなる.その結果,結晶弾性の成分が生じてくる大変形領域では少しのひずみの増大で応力が急激に増大する.このように考えると,粘着剤の弾50

性率がひずみの関数となり,ひずみの増大によって指数関数的に弾性率が高くなる
ことも容易にうなずける.

4.3.2 速度形定式化
ここでは有限要素法における粘着剤の数値シミュレーションのために,上記で述
べた材料モデルの3次元状態での速度形の構成方程式を導出する.いま3次元状態
でのひずみEのノルムe‑を次式のように定義する.

EA ‑斤寺(4.40)

そして3次元状態における式(4.39)のEiを, 6‑を用いて次式のように設定する.
E% ‑ Al exp(Ble)(4.41)

いまi番目要素におけるスプリング部分のひずみをE皇pとすると, i番目要素におけるCauchy応力qiの偏差成分g′『こついての構成方程式は次式のようになる.
g′

‑

2Gle芸

(4.42)

本来は式(4.42)の右辺のひずみについても偏差成分を用いるべきであるが,本モデ
ルでのスプリング部分は非圧縮と見なせるため,通常のひずみを用いている.ここで, GiとEiとの関係式を再掲しておく.
E%

Gr‑(4.43)

2(1 +1/*式(4.42)の物質微分を考えると次式のようになる.

‑ 2G莞里+2芸'sp
DGl(4.44)o‑′i= 2GID呈p + ‑石『存ここでD呈pはi番目要素におけるスプリングのひずみ速度である・またここでは式(4.44)のGiも時間の関数であるので,その物質微分は式(4.41)(4.43)より次のようになる.
DGI AI D r
Dt

2(1+vl‑ BIGl芸51

義{exp郎)(4.45)

構成方程式を大変形問題にも対応できるようにするために, q′iの物質微分は客観応
力速度へ置き換えられるものとする囲.ここではいくつかの客観応力速度の中でも, Jaumann速度を用いる.その場合のJ′iのJaumann速度は,式(4.44)(4.45)を考慮して以下のようになる.

(筈(J) ‑ 2(PDiBp+Bi^(4.46)

そしてダッシュポット部分の構成方程式は一般化Maxwellモデルの時と同様に次式となる.
crn ‑ 2り蝣D¥(4.47)

Fig.4.2のモデルのひずみ速度Dは,スプリング部分Dspとダッシュポット部分Ddpの和で表現できや.さらに上記でも述べたようにモデルのすべての要素における速度が同一であるので, Dもすべての要素で同一である.よってモデル全体でのひずみ速度βが次式のようになる.
D ‑ D呈p+D昌(4.<

式(4.46)(4.47)(4.48)を整理すると, i番目要素におけるq′iの速度形構成方程式は次のようになる.

(笠(J) ‑2G{D+B^a*‑げl2(1 + z/jr*(4.49)

cauchy応力の偏差成分q/iの速度形構成方程式が材料モデルの要素ごとに定義されるのに対して,圧力pに関してはすべての要素を含めた線形弾性体として定義する.その場合の圧力pは次式のようになる.
p ‑ ‑Kvtie(4.50)

ここで体積弾性係数Kvは次のようになる.

Kv‑差〈

El
(4.51)

一般化Maxwellモデルの場合とは異なり,修正モデルでは体積弾性率Kvは時間の関数となる.よって,式(4.50)の物質微分を考えると,一等tre ‑ KvtvD

墜p‑K,,trD
Dt Kv
52(4.52)

となる. Kvの物質微分を考えると,
B%E%
1.53)

となる.式(4.52)(4.53)より圧力pの物質微分は以下のように書き換えられる.
‑ K,,trD

禁‑孟差〈(4.54)

圧力pはスカラーでありスピンWの影響を受けないので,式(4.54)は大変形状態においても成立する.
以上をまとめると, Fig.4.2に示すような本研究で修正した一般化Maxwellモデルの大変形状態にも対応した3次元状態での速度形構成方程式は

a‑&‑pZを考慮して次式のようになる.

2GID+B瀞i ‑ql・2(1 +1/*: ・・.I3(1 ‑ 2i/')

53

芸‑KvtrD I (4.55)

第5章

3次元PIJIC法による界面捕捉

粘着剤は容易に大変形を生じるため,その変形を記述する解法としてEuler型解法を用いている. Euler型解法は大変形の記述には有利であるが,物質境界面の捕捉精度と境界条件の取り扱いに課題を残している.通常, Euler型解法では界面の捕捉には3.5章で述べたVOF法が幅広く用いられてきたが,長時間積分においては物質境
界面の数値拡散により界面捕捉の精度が低下してしまう.また, Euler型解法では,Lagrange型解法のように物質境界面が明確に定義されておらず,境界面上に格子の節点がないため,境界条件を与えることが困難である.つまり,一般にLagrange型
解法で行われているように,物質境界面に対して力を与えたり(力学的境界条件) ,強制変位を与えたり(幾何学的境界条件)する問題を取り扱うことができない.このことによりシミュレーションで記述できる変形の形態が限定している. VOF法による界面捕捉法では, ALE法を用いることによりかろうじて‑軸引張り解析を行う
ことはできるが,ピールなどの実用的な解析を行うことができない.
そこで, 5章と6章では任意の境界条件を付与する手法について説明する.まず本章で, Euler型解法にて物質境界面を定義する手法について説明する.次章において,本章で説明した手法を用いて定義された境界面に境界条件を付与する手法について説明する.
Euler型解法にて物質境界面を定義する手法として, 3次元PLIC法による界面捕捉法を用いている. 3次元PLIC法は. (1)界面の法線ベクトルの計算. (2)境界面の再構築. (3)境界面の移流の3つ計算手順に分けることができる.法線ベクトルの計
質では,注目する界面セルの周囲の体積率分布から計算する比較的簡単な方法を用いている.境界面の再構築では,従来のPLIC法で用いられることの多かった反復計算を必要とする計算コストの大きい手法でなく,より高速で高精度な直接計算法を用いている.境界面の移流では,アルゴリズムが簡便になる各座標軸方向の移流を別々に行う方法を用いている.

54

5.1

3次元PIJIC法の概要

PLIC法は,各セル内で物質境界面を線形近似してセル境界を横切る立体の体積を計算し, VOF関数を移流するL=とによって,移動境界面を高精度に捕捉する手法である. HirtとNichols岡により提案されたVOF法(1981)は,移動する物質境界に対する数値解析手法として,これまで幅広く用いられてきた.近年はVOF法に代わって,より少ない格子分割数で高精度な結果が得られるPLIC法の適用が増えつつある. PLIC法は, Youngs[8]による先駆的な2次元PLIC法(1982)に始まり,さらに高精度な2次元PLIC法は1980年代から1990年代にかけて数多く開発されてきた[76][77].一方, 3次元PLIC法は,幾何学的な演算の複雑ざから開発・適用事例は多くはなかったが, 2000年代に入って,高精度な3次元PLIC法の開発・適用事例が増えているlgト[15].
本研究の3次元PLIC法のアルゴリズムは, Youngsが提案した2次元PLIC法(詳細は陶に詳しい)及びScardovelliとZaleski閏の提案した計算法に基づいている.
PLIC法では,線形近似した物質境界面の方程式を決定するために,平面の法線ベクトルと原点から平面までの距離を求める必要がある.法線ベクトルの計算法には,Chorin[80], Barth[81], Swartz[82]などによる方法が挙げられる.本研究では,計算コストの小さいYoungsの計算法岡を採用している.また,原点から平面までの距離は, Kotheら[76]のように通常は距離に関する3次方程式の解をBrent法国などの反復計算により求める.一万,本研究ではScardovelliとZaleski閏が提案した, 3次方程式の解が陽に表された代数式により計算を行うため,反復計算は不要となり計算コストが大幅に軽減される.また,移流項の差分では,アルゴリズムが簡便なdirection‑split法を用いている.このように,本研究の3次元PLIC法のアル
ゴリズム特徴は,計算コストの少なさとdirection‑split法によるアルゴリズムの簡便さにある.

55

5.2 境界面の再構成
5.2.1 境界面セルの判定
xl,>x2X¥

Fig.5.1: 0<￠<1

Fig. 5.3: ￠‑0

Fig.5.2: ￠‑1

Fig.5.1のように,物質の境界面が存在するセルを,境界面セルと呼ぶことにする.
セルにおける物質の体積率を￠としたとき

0<6<1を満たすセルを境界面ゼルと判定し,境界面セルのみで境界面を再構成する.実装上は,
E≦4,≦1‑E(5.1)

を満たすセルを境界面セルと判定する. Eには e‑1.0×10 3の値を与えている.
体積率1のセル(Fig.5.2)と体積率Oのセル(Fig.5.3)の場合は,一次風上差分法により体積フラックスの移流計算を行うため,境界面の再構成は不要である.

5.2.2 線形近似された境界面の方程式
直交デカルト座標系(xi,x2,x3)において, Fig.5.4に示す辺長△3:iの直交格子を考える.線形近似された境界面の法線ベクトルをn‑ (ni,n2,n3)とすると,境界面の方程式はriiXi +112X2 +n3二r3 ‑ α(5.2)

で与えられる.ここで, αは定数であり,原点から平面までの距離と関係付けられる1.式(5.2)より,境界面を再構成するには,各セルにおいて線形近似された境界1原点から平面までの距離をdとしたとき,点と平面の距離の公式より, dは次式で与えられる.
d=
＼ ∴

/..

56

/蝣‑.

Fig. 5.4: Linear interface in each cell面の法線ベクトルnと,距離定数αを求める必要があることがわかる.
ここで,次式のような関数を定義しておく.
f(xl,x2,x3) ‑nlxl+n2x2+n3x3‑α(5‑3)

この関数f(xl,x2,x3)は,境界面の法線ベクトル方向の領域では月xi,x2,x3) > 0となり,一方,逆向きの領域ではf(xl,x2,x3)<0となる.

5.2.3 境界面の法線ベクトルの計算
境界面の法線ベクトルnは,外向き方向を正にとった場合,各セルの体積率(VOF関数) ￠より以下の式で与えられる.
n‑‑∇￠

(5.4)

式(5.4)の体積率勾配の計算法として,本研究ではYoungsの方法閏を用いる.他の計算法として, Chorin[80] (l次精度, 2次元のみ), Barth[81] (l次精度), Swartz[82](2次精度)によって提案・された方法が挙げられる.一方, Youngsの方法は1次精度で,これらの手法に比べ3倍から5倍ほど低い計算コストで済む手法である[叫.

57Youngsの方法では(ij,k)の位置のセルにおける体積率勾配を以下のように差分近似する.

(芸i,j,k4>i+l,j,k ‑ 4>i‑l,j,k￠i,j+l,k ‑ Oi,j‑l,fe2△￡2dxi i,j,k・fii,j,k+l ‑ Oi,j,fc‑1

2△3:3

dx? i,j,k

ここで￠は,次式で定義されるように,(5.5)

JA.′・l

(5.6)

5.7)

セル%hk)隣接するセルの体積率をパラ

メータa,b,cにより重みつき平均した体積率である.

串i+l,j,k‑[a{十l,j‑l,k‑1+Oi+l,j‑l,fe+l+4>i+l,j+l,k‑1+Oi+l,j+l,fc+l)+K(f)i+l,j‑hk+4>i+l,j+l,k+(Pi+l,j,k‑1+^i+l,j,fe+l,+C(pi+1Jtk]/(4:a+4b+c)(5.8)

串i‑l,j,k‑[a(れl,j‑l,k‑1+叫7‑1,/c+l+(Pi‑lj+l^‑1+<^i‑lj+l,fe+l)+K4叫j‑l,k+￠叫j+l,k+<Pi‑l,j,k‑1+￠i‑lj,fc+l,‑l,3,k}/(4a+4b+c)(5.9)∂i,j+l,k‑[&(</>i‑lJ+l,fc‑1+^i‑lj+l,fc+l+Oi+l,j+l,/c‑1+(Pi+lj+lfi+l,+b((/)i‑itj+itk+￠汁l,j+l,fe+Oij+l,fe‑1+Oi,j+l,fc+l,+cOiJ+iifc]/(4a+46+c)(5.10)∂i,j‑l,k‑[&(<&‑!,上1,/c‑1+<!>i‑l,j‑l,k+l+4>i+l,j‑l,k‑1+￠i+lj‑l,/H‑l,+6(Oi‑ij‑i)fe+cj)i+ij‑itk+Oi,j‑i,fc‑i+</>i,j‑i,fc+i,+cOij‑‑i,fc]/(4a+46+c)(5.ll)

示,j,k+l‑[a{4>i‑l,j‑l,fc+l￠i‑l,j+l,fc+l+4>i+i,j‑l,fc+l+Oi+l,j+l,fc+l)+b(￠i‑l,j,k+l+Oi+lj,fc+l+￠i,i‑l,fc+l+<Pi,j+l,k+l)+C<t>i,j,,.fc+lj/(4a+46+c)(5.12)∂i,j,k‑l‑[a{(Pi‑l,j‑l,k‑1+￠蝣i‑l,j+l,k‑1+(j>i+l;j‑l,k‑1+<f>i+lj+l,k‑1)

58

+b(￠i‑l,j,k‑1+4>i+l,j,k‑1+4>i,j‑l,k‑1+</>ij+l,fc‑l,+C(f>ix%j,k‑1]/(4a+46+c)(5.13)式(5.8)(5.9)(5.10)(5.ll)(5.12)(5.13)からわかるように,Youngsの方法では,セルh3,k)の勾配を単なる中央差分法により計算するのではなく,セルi,j,k)に隣接す
るセルにステンシルを広く取り,重みつきの平均値で差分近似を行っている.Youngs
は,パラメータa,6,cについて次の値を使用している.

a‑2,

6‑2,

c‑4

(5.14)

6‑0,

c‑l(5.15)

ところで,各パラメータの値がa‑0,の場合は,セル'ij,k)における体積率勾配は(芸i,j,k(芸i,j,kW+lj.fc ‑ <Pi‑l,j,/c2△￡1

vi,j+l,k ‑ <Pi,j‑l,k(5.16)(5.17)

・̲'A./蝣

4>i,j,k+l ‑ ￠i,j,k‑1

2△3;3・dx* i,j,k(5.18)となり,中央差分法に一致することがわかる.

5.2.4 距離定数の計算
次に,距離定数αの計算法を説明する.いま,境界面の法線ベクトルnの成分が全て正である場合のみを考える.このとき, Fig.5.5に示すような辺長△3:iの直方体(セル)によって切り取られる立体(ABGH‑LMNK)の体積のま,以下で与えられる.
3

6nin2n3‑∑F3(α‑ni△zi)+∑F3(α‑α責F3(a‑α +714△zi)+riiAxi) (5.19)

α31蔓

i=l

式(5.19)において, αmaxと関数Fnをそれぞれ次式のように定義している.
3・max ‑ ∑ni△Ti
1‑159(5.20)

Fig. 5.5: The truncation volume ABGH‑LMNK

for y>OFn(y)‑ y"(5.21)for y≦0

式(5.19)の第一項α3/6nin2n3は四面体AUKの体積である.第二項∑=1*3 α‑ni△Xi /6ni712^3は,頂点I,J,Kがセル境界を越えた場合,四面体AUKから取り除かれる四面体の体積である. Fig.5.5の場合であれば,四面体HIPL, BOJNである.第三項∑Ll^(α‑α+ui△zi)/6nin2n3まFig.5.5の四面体GOPMのよ

うに,直線IJ, JK, KIがセル内を通過しない場合に第二項が重複して取り除く
四面体の体積である.
ここで,セルの辺長とαmaxを正規化した場合を考える.セルの辺長については,△￡」 ‑ 1

(5.22)の場合を考える.また, αmaxについては,平面の方程式(5.2)の両辺を∑1=1Hiで割ることにより,
3・max‑ ∑us ‑ 1(5.23)

i‑lと正規化する.以上の場合,体積Ⅴ及び距離定数αは以下のような性質を持つ.
60

(1) Vつま距離定数αのCl級関数で,単調に増加する.
(2) Vとαの値の範囲は[0,1]である.
(3)のまnun2,n3の順列に対して不変である.
(4) Vとαのグラフ(一例をFig.5.6に示す)は,点(V,α ‑(1/2,1/2)に関して点対称である.

e
e
'J=
U
Cq

と
Q)

≡
s
plllllO>

Fig. 5.6: The graph of Volume fraction

性質(2)より,体積Vと体積率(VOF関数)￠は一致することがわかる.また,性質(3)よりniの順列について以下の場合のみ考えることにする.
0<71i≦n2≦n3(5.24)

さらに,性質(4)より, Vとαの値のおいて次式の範囲のみを考えればよい.
0<V<(5.25)
10<α<‑2(5.26)

ここで,式(5.22)(5.23) (5.24)が成り立つとき, αを陽な式で表現するために, αの値の範囲により式(5.19)を場合分けして簡略化する.式(5.24)が成り立つ場合,
61

Fig.5.7: 0<nx≦n2≦n3四面体AUKの各辺長はAK<AJ<AIとなることから,四面体AUKはFig.5.7に示すような立体になる.よって, αの値によって後述するように場合分けすることができる.

62

(1)0≦α<nl

Fig.5.8: 0≦α≦nl

四面体AUKがFig.5.8のような場合,頂点月まセル内に位置する.したがって,この場合のαの制約は次式になる.

o<^<l
mil

⇔0

≦α

<nl

(5.27)

式(5.27)が成り立つ場合,式(5.19)より体積Vは次式になる.
V=・I

6nin2n363(5.28)

(2)m≦α<n2Fig.5.9: nl≦α<n2四面体AUKがFig.5.9のような場合,頂点1はセル外に位置し,頂点J, Kはセル内に位置する.したがって,この場合のαの制約は次式になる.

⇔nl ≦

α

<n2

(5.29)

式(5.29)が成り立つ場合,式(5.19)より体積Vは次式になる.

6nin2n3[α3 ‑ (α‑no3](5.30)α(α ‑ n¥)+n至2n2ri36n2n364

(3) n2≦α<nFig.5.10: n2≦α<n四面体AUKがFig.5.10のような場合,頂点1, JLまセル外に,頂点Kはセル内に位置し,直線IJはセル内を通過している.したがって,この場合のαの制約は次
式になる.

α
‑

<

n3

it‑i+n2‑α>0

⇔n2

≦

α

<n

(5.31)

ここでnは, nl+n2とn3のうち小さい方の値であり,次式により定義する.
ォi2

n

‑

‑

ni+no

mm(n12,n3]

式(5.31)が成り立つ場合,式(5.19)より体積1月ま次式になる.

6nin2n3[α3‑(α‑nx)d‑ α‑n2f]

α!(3n12‑α +nf(ni‑3α)+n召(n2‑3α)
6nin2n365(5.34)

(4) n≦α≦1/2n≦α≦1/2の場合は,以下の2通りの場合が考えられる.
(4‑1) n3≦α≦1/2

Fig. 5.ll: n3≦α≦1/2n3≦α≦1/2が成り立つ場合,四面体AIJ狛まFig.5.11のように,頂射, J, K
はセル外に位置する.また,直線IJuセル内を通過する.なぜなら,この場合のαの制約は次式になるためである.
α

1
n3

n3三

ni+n2
1

⇔n3≦α≦亘

ここで(5.35)

n¥+n2の値の範囲について考える. n3≦ n¥+n2とn¥+n2+n3‑1が成り立つことにより,
1‑(ni+n2)≦ni+n2

芸≦ n¥ +n2となる.いま, 0≦α≦1/2の範囲を考えているため(5.36)
n¥+n2≦αについては考えなくてよいことになる.よって,式(5.35)が成り立つ場合,式(5.19)より体積のま66

次式になる.
V=
6nin2n3[α3‑(α‑niY‑(α‑n2y

α‑nsf]

α!(3‑2α)+n至(ni‑3α +nl(n2‑3α +n￨(n3‑ α)(5.37)

6r乙i n2n3(4‑2) n12≦α≦1/2

Fig. 5.12: nu≦α≦1/2n12≦α≦1/2が成り立っ場合,四面体AUKがFig.5.12のように,頂点1, Jはセル外に,頂点Kはセル内に位置する.さらに,直線IJはセル外を通過する.なぜなら,この場合のαの制約は次式になるためである.
n¥+n2‑α≦O(rii+n2 ≦n3
1⇔n12≦α≦亘(5.38)

ここで, n3の値の範囲について考える. nl+n2≦n3とnl+n2+n3‑1が成り立つことにより,
1‑n3<n31≡‑":

67(5.39)

となる.いま, o≦α≦1/2の範囲を考えているため, n3≦αについては考えなくてよいことになる.したがって,式(5.38)が成り立つ場合,式(5.19)より体積Vは次式になる.
6nin2n3[α3‑(α‑ni) ‑ α‑n2f+ α‑nx‑n2f](5.40)

2α‑n122n,以上の体積Vについての場合分け(1). (2). (3), (4)を以下にまとめる.
v‑

α3

forO<α < nl
6nin2n3v=

I//=α(α ‑ rii) n至+
2n2U3
Qn2fi3

forni≦α < n2

α!(3ni2‑α)+n亨(nx‑3α +rq(n2‑3α)oni^r^

for712≦α < nα2(3‑2α +n?(ni‑3α)+n∃(n2‑3α +nl(n3‑α)

forn^≦α≦ 1/2(5.44)

6nin2n32α‑n12for ni2≦α≦1/2n,(5.45)

以上の式(5.41)‑(5.45)を変形することより, αは以下のように表せる.ただし,
0≦α≦1/2であることを考慮している.

6nininsV for O≦V<Vl
n至+8n2n3(^‑ ^;(5.46)
for Vx ≦ V< V2(5.47)P(α‑<2qα +a'9α +a[α+a乙‑0 for V2≦v<v*(5.48)P(α‑<2。α3+a"α2+a'lα+an‑O for V331≦V≦1/(5.49)

・ ‑ n3V+警for V32≦V≦1/(5.50)

式(5.46)‑(5.50)において, vlf v2, v31, v32, v3を以下のように定義する.
Vl‑V ¥ 。=niV^v¥α=花‑max(6n2ri3, e)
Vi

+

(n2

‑

nl)/2n3vzl‑V¥。=れ‑ [nl(3n12 ‑n3) +n冒(ni ‑ 3n3) +n2(n2 ‑ 3n3)] /6nin2n3 (5.53)

Vz2‑Vl。=ni2 ‑ n12/2n3
vs‑v

α=n

‑

mm(V31,(5.54)

V32)(5.55)bS

ここでVlは ni ‑n‑2 ‑0の場合の零割りを避けるためのn冨/6712^3の実装上の近似式である. Eには零でない,十分に小さな値を与える.さらに, αに関する3次方程式(5.48)(5.49)の係数はそれぞれ以下のようになる.
り:;
/

I

(5.56)

3n12(5.57)‑3(n至+n芸)(5.58)

n号+nヨーQninznzV(5.59)

‑

‑

a[ ‑
/

‑

//

1/. ‑

ご

(5.60)<2r, ‑

3(5.61)a'[ ‑‑3(n至+ni+n…)(5.62)

n至+n芸+ni‑6nin2ri3V(5.63)//‑

αに関する3次方程式(5.48)(5.49)には, 3つの実数解が存在し,求めるべき解は2番目に大きい解になる囲.この3次方程式の方程式の厳密解は,以下のようにして求められる囲.まず, α3‑ 1とするために3次方程式(5.48)(5.49)の両辺をα3で割り,新たにpo> <?Oを以下のように定義する.
α1 α∃
3

9

a¥a,2‑3ao a塁
627

さらに, βを次式で定義する.
∂=王C。S‑1
3(5.66)

以上により,求めるべき解は次式で与えられる.
α ‑ √二元(諦sin6 ‑ cos#ト聖3(5.67)

5.2.5 境界面の場合分け
前節では,境界面の法線ベクトルnの各成分が次式を満たす場合のαの計算法を説明した.
0≦nl≦n2≦n3
69(5.68)

その他の場合は,式(5.68)の場合に変換してαを求める.まず, Fig.5.13に示す
ように, niの正負により8通りの境界面が考えられる.さらに Fig.5.14に示すよ
うに, niの絶対値の順列で6通りの境界面がある.つまり,法線ベクトルの成分の値により,合計48通りの境界面が考えられる.

70

Ai

l m≧0,n2≧0,n3≧(2)nx≦0,n2≧0,n3≧O

X3

X3

Xl
(3)ni ≦0,n2≦0,n3≧(4)711≧0,n2≦0,713≧O

X3

X3

A"i

(5)ni≧0,n2≧0,n3≦

6 ni≦0,n2≧0,n3≦O

X3

X3

Xl
(7)m≦0,n2≦0,n3≦(8)m≧0,n2≦0,n3≦O

Fig. 5.13: Interfaces in case n; is positive / negative71

(l)m≦Kl≦H(2)nx ≦Kl≦H
X3

X3(4)￨n2 ≦*3 ≦n¥¥(3)n2≦ni ≦n3

A3

X3(6)恒3l< n2 ≦hi(5)¥n3¥≦kl≦H

Fig. 5.14: Interfaces by the order of the absolute value of rij

72

5.3 境界面の移流
5.3.1 移流方程式の離散化
前節までで,境界面の再構築について説明した.本節では,再構築した境界面の移流計算の方法について説明する.体積率(VOF関数) ￠の移流方程式は,

霊+(v ∇)桓
で与えられる.連続の式∇(5.69)

v‑Oが成り立つことにより,発散型に変形する.

霊+,・(V￠)‑o(5.70)

ここで,体積率のフラックス′を/‑ォ<;(5.71)

霊+芸+芸+芸‑0(5.72)

とおくと,移流方程式はと表せる.式(5.72)を時間方向に前進差分近似すると・Jl+l ￠和一△砦‑△萱‑△萱

(5.73)となる.さらに空間方向に差分近似することにより,
LU+1‑鶴一芸(恥,k萄,3,k.芸(fn‑」n¥‑芸te¥<">jifc+i・n
i,j,k一芸)(5.74)となる.さらに,境界面が1ステップでセル境界を横切る体積(体積フラックス)を求めて,移流計算を行うために,式(5.74)を次式にように変形する.

媚+1‑jl)n.‑△yn
i+吉,j了△△i,3+去,k‑△yn去,kVii,j,kVii,j,k△yn
i,j,k+去‑△吉Vii,j,k(5.75)

式(5.75)において Vihj,kはセルの体積であり,次式で定義される.
Vii,j,k △∬1△￡2△3;373(5.76)

また,体積フラックス△Vはそれぞれ以下のように定義される.
△yn
i+去,j,k瑞,j,k:△z2△￡3△t(5.77)△vn‑吉J>k‑tnuk△T2△T3△t(5.78)△yn吉‑fn去,k△zl△z3△t(5.79)△yn寸m寸k△zl△￡3△t(5.80)△v:rnij,fc+￨‑‑<>i,j,k+吉△El△￡2△t(5.81)△去‑fn.
ijfe一言△zl△z2△t(5.82)

さらに,アルゴリズムの簡便さを考慮して,以下のように式(5.75)を方向分離して計算する.

鶴k‑鶴ki,j,k bh,k‑

kn+1 ‑砿k△yn
i+圭,j,k△(5.83)

Vii,j,k△V*
1,3+吉,k‑△V*.
*,3‑去,kVii,j, k(5.84)‑11∴ ‑Al"(5.85)

Vii,j, k

5.3.2 体積フラックスの計算
はじめに全てのセルにおいて1次風上差分法により体積フラックスを計算する.
次に,境界面セルにおいて式(5.41)‑(5.45)を利用して体積フラックスを再計算する.その際,セルから流出する方向の体積フラックスの場合のみ, 1次風上差分法
により計算した体積フラックスを更新する.

74第6章 境界条件の付与と空間離散化
本章では, PLIC法により得られた物質境界面に対して境界条件を付与する方揺について説明する.境界条件の付与は基礎方程式の空間離散化の際に行う.空間離散化には,有限要素法と差分法を用いている.有限要素法を用いることにより,力学的境界条件及び幾何学的境界条件を取り扱いやすくなる.本研究では幾何学的境界条件の定式化にはペナルティ法を用いている.

6.1 境界条件の付与
本節では,境界条件の付与法について説明する. 3次元ユークリッド空間R3内において,ある物質が占める領域n⊂R3を考える. Fig.6.1に示すように,領域口の境界Il‑∂0の一部分Ilf ⊂ Ilにおいて表面力が与えられているものとする.また,境界の一部rlu ⊂ rにおいて幾何学的境界条件が課せられているものとする.

Fig. 6.1: The domain and the boundary of moving body

このとき, Lagrange表示の平衡方程式(3.22)を再掲する.
paL

∇

75

cr+pb

(6.1)

また,力学的境界条件,幾何学的境界条件はそれぞれ以下のように表される.
cr・n‑tonr,u‑u on F,:

ここで, nは境界上で定義される外向き単位法線ベクトルである.また,占t, uはそれぞれ,既知の情報として与えられる物体力ベクトル,表面カベクトル,変位
拘束ベクトルである.
ただし,領域境界IlfとIluは互いに重なりを持たず,
rtnru‑功かつ

r‑r.Uru(6.4)

が成り立つことを仮定する.なお,∩は積集合,鋸ま空集合,∪は和集合を表す数学記号である.
平衡方程式の両辺に任意の仮想変位6uを内積し,領域0で積分した式(3.27)を再掲する。
/paL蝣SudV+cr:(5u'aJn㊨∇)dV‑fJn∇・‑5u)dV+pb‑5udV(6.5)
式(3.27)から仮想仕事式(3.28)を導く時にGaussの発散定理を用いるが,そのJq
際,力学的境界条件(6.2)を用いることにより,/o∇su)dV‑I'n‑{cr‑5u)dS
Jy

(6.6)
t蝣5udS

Ttとなる.式(6.6)を式(6.5)に代入することにより,以下の式を得る.
pa,L蝣SudV+cr:(Su㊥∇)dV‑f
Jv予∂udS+/pb‑SudV(6.7)
Ju
これは仮想仕事式であるが,力学的境界条件(6.2)は,式(6.6)に示すように,仮想仕事式(6.7)の導出の過程で取り込まれている.一方,式(6.3)で表される幾何学
的境界条件は,仮想仕事式(6.7)に取り込まれない.したがって,仮想仕事式を有限要素法により離散化した方程式を解く際に,幾何学的境界条件を付帯条件として与える必要がある.そこで,本研究では処理が単純なペナルティ法[86]を用いて幾何学的境界条件を処理することにする.

76

ペナルティ法により幾何学的境界条件を付加した仮想仕事式は以下のようになる.

paL蝣SudV+ / a: (5u<g>∇)dV+ITu戸(u‑由蝣5udS‑ / t‑SudS+ pb‑SudV(6.8)

ここで,密はペナルティ係数である.通常は固体の弾性係数の104倍程度の実数を設定することが多い.

6.1.1 有限要素法による空間離散化
ペナルティ項を付加した仮想仕事式(6.8)に対して,有限要素法により空間離散化を施す.いま,領域0をNel個の小領域Oeに分割することを考える.この小領域0。
は有限要素と呼ばれる.このとき,各要素領域Oeは次の関係を満たすものとする.
JVe1

0‑Uoe
c

かつ

n。1nne2‑め(el≠ea)

(6.9)

1

式(6.9)に基づく分割は,領域0を隙間や重なりがなく,要素で埋め尽くされるよ
うな分割を意味する.各有限要素領域0。においても仮想仕事式(6.8)は成り立っため,仮想仕事式(6.8)を以下のように表すことができる.

射pcLl・SudVe善上a:(su ∇)dVe+ Nei」e=l''Tue戸(u‑面SudSp
iVel

岩上e

t蝣8udSe

e=l''Tiepb蝣6udVee(6.10)

式(6.10)における加速度aL,仮想変位6u,変位uをそれぞれ,各要素の節点値&l, 5u, uにより以下のように近似する.
Ill

̲ヽ "/

6u竺NdutL.三Nl‑I.

ここで, Nは重み関数から構成される形状関数マトリクスである.本研究では,有限要素としてFig.3.2に示す3次元8節点アイソパラメトリック要素を用いる阿国.

77

次に,各要素領域の和で表された仮想仕事式(6.10)に,式(6.ll)(6.12)(6.13)を代入することにより,以下の式を得る.
MaL + Fint + Fp ‑ FEext
ここで, Mは質量マトリクス(6.14)

Fintは内カマトリクス, Fpはペナルティ法による付加拘束マトリクス ・extは外力マトリクスであり,それぞれ以下の式で与えられる.

M=

‑Tint

NTNdVr(6.15)¥l Bl crdVe(6.16)

善β上e‑‑?E/r(nTn(6.17)
e=lJTu崩‑NTu)dSe
Frext善上NTidSe+J2pNTbdVe(6.18)
te=iJneここで, Nuは幾何学的境界条件が課される有限要素の数, Ntは力学的境界条件が課される有限要素の数である.式(6.18)において,右辺第一項は表面カマトリクス
であり,この項によって力学的境界条件が付加される.また,式(6.16)においてBLは変位‑ひずみ変換マトリクスであり,以下のように定義される.
0

0

0

o o

慧

0

0

0

o o

0

∂Ⅳ1

∂N2

∂N8

∂y

ay

∂y

O

O

O

・JVt

∂N2

∂N8

手て

dz

′1.

6.1.2 選択型低減積分
本研究では,非圧縮性や微圧縮性の物質を取り扱うため,体積一定の拘束を課していることになる.有限要素解析では,この体積一定の拘束を各要素の積分点ごと78

に評価するため,変形が本来よりも拘束されて結果が剛となりすぎる体積ロッキン
グ現象が起こることがある.この体積ロッキング現象を避けるために,本手法では選択型低減積分法を用いる.この方法は速度勾配テンソルを偏差成分と体積成分とに分割して, Fig.6.2のように,体積ロッキングに関連している体積成分のみを低減
積分し,残りの偏差成分については完全積分を施すというものである.

△ Integration points for volumetric strain
ロIntegration points for deviatoric strain

Fig.

6.2:

Integraノtion

Pointsまず,以下のように速度勾配テンソルLを偏差成分L′と体積成分Lvに分割する.
L=L'+Lv(6.20)

ここで,体積成分は以下のようである.

・V 吉tri/(6.21)

いま,低減積分において評価した速度勾配テンソルを'redとする.選択型低減積分
では,体積成分だけ低減積分を行い,残りの偏差成分は完全積分を行うわけであるから,選択型低減積分における速度勾配テンソルLselは完全積分における速度勾配テンソルをLfulとして以下のように評価できる.

Lsel = Lful + I ^tlLred I‑ (hrLful(6.22)

上式は体積ひずみに関与する垂直ひずみ成分のみ低減積分を行い,ひずみの偏差成
分は完全積分することを意味している.これを離散化した有限要素解析で反映させ79

るためには,以下のようにBLをBfのように修正すればよい.
(6.23)

Bf‑B[ul+吉Tired芸蝣‑ful
ここで,

BLv ‑

∂x dy

∂Z
0

∂x dy
0

∂2:

0

0

0

0

0

0

0

0

0

0

...

0

0

0

0

0

0

0

0

0

...

0

0

0

である.上式のEjsetを用いて式(6.16)の内カマトリクスを求める.
構成方程式で用いる新たなマトリクスBNLについてもここで定義しておく

BNLを以下のように与える.
0

0

0

0

0

0

0

0

0

0

∂z
0

BNL ‑

0

∂2:

0

0

0

0

0

0

0

0

0

0

0

0

∂z

∂z

∂Nl

0

∂z

dNo

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0
∂z

∂Ⅳ1

0

0

このBNLについても, BLと同じように修正してf>selを算出する.すなわち

TDsel

を以下のように求める.
‑tsel‑>ful.吉B鑑一芸ptful
‑NLv

80

(6.26)

ここで,!ノ

BnLv ‑

<hi

,I二

0

0

0

0

0

0

0

0

0

∂N8

∂N8

∂N8

l).r

ay

∂z

O

O

O

O

O

O

0

0

0

0

O

O

O

∂N2

dN*

∂N8

∂N8

∂N2
・hi

,)‑

否㌃

dz

(6.27)と与えられる.以後,ひずみの評価や応力の評価で用いる全てのBL,BLの代わりにB芸ITDsel
NLを用いることで,計算の整合性を持たす.

6.1.3 力学的境界条件の処理
ここでは,力学的境界条件を付加する表面カマトリクスについて説明する.いま,外力マトリクス(6.18)を次式のように表すことにする.
Fenct.

Ft(6.28)

ここで, Ftは表面カマトリクス, Fbは物体カマトリクスであり,それぞれ次式で表す.
Nl tdSp

*t‑

善上土

Fh‑

至上e pNl bdVe

(6.29)(6.30)

表面カマトリクス(6.29)について,各有限要素において面積分を行うのは計算コストが大きい.よって本研究では, Fig.6.3に示すように,要素内における表面力が,PLIC法により定義される境界面の重心に集中荷重として作用しているものとして取り扱う.重心の算出の際には, PLIC法により求められる各要素の境界面の方程式81

Fig. 6.3: Traction at piecewise linear interface

から,要素の各辺と境界面との交点の座標を初めに計算する.次にこの交点を頂点とする多角形の重心を計算することで境界面の重心が決定される.このように表面力を重心に作用する集中荷重として取り扱うことにより,数値積分を行う必要がな
くなり,簡便な処理で済む.この場合,表面カマトリクスは次式で表せる.
=＼11

Ft ‑ ∑NTtSp(6.31)

蝣=]式(6.31)において, Seは各要素における線形近似された境界面の面積である.式(6.31)を成分表示すれば,以下のように表せる.
iVi
0

0

0

Ni

0

0

Ⅳ1

m
ty

N8

0

0

0

JV8

0

0

0

iVs

sp

tz

(6.32)

Nitx

Wy
iVt

Nxtz

≡
N8 tx
Nsh
AW8 <<z82

式(6.32)における各節点の重み関数Mi ‑ 1‑8)は,局所座標系における境界面の重心の座標vsg>vgiS3を代入することにより,次式で与えられる.

・1‑芸 (1‑&)(1‑%)(1‑Eg)
N,=‑ (1+Q(l‑vg)(l‑Q
N*=‑ (1+Q(l+vg)(l‑Cg
NA‑‑ (1‑&)(i+サfc)(i‑Cサ(6.33)

・5‑芸(1‑Q(l‑vg)(l+Cg

・6‑芸(I+Q(l‑vg)(l+(g・7‑芸(1+&)(1+%)(1+Cff
・8‑芸(i‑&)(i+サfc)(i+C,)

6.1.4 幾何学的境界条件の処理

Fig. 6.4: Displacement at piecewise linear interface

ここでは,幾何学的境界条件を付加するペナルティ法による付加拘束マトリクス(6.17)について説明する.ペナルティ法による付加拘束マトリクス(6.17)について
各有限要素において面積分を行うのは計算コストが大きい.よって,力学的境界条
件の処理と同様にして, Fig.6.4に示すように各要素の境界面の重心に幾何学的境界条件が課せられているものとして取り扱う.この場合,ペナルティ法による付加拘束マトリクス(6.17)は次式のように表せる.

83

Wfl

Fp ‑p∑ (NTNu ‑ NT面) se(6.34)

・=1

Euler型解法においては,メッシュは空間に固定されており,節点変位は計算されないため,式(6.34)において節点変位及び変位拘束ベクトルは,それぞれ以下のよ
うに計算する.
U=V△t
u=V△t

ここで, △自ま時間増分, i)は前ステップ節点速度, i)は既知の情報として与える速度拘束ベクトルである.式(6.35)(6.36)をペナルティ法による付加拘束マトリクス(6.34)に代入し成分表示すれば,以下のようになる.

ffiff

N,

0

0

0

iVi 0

0

0

JVi

Fp‑pE
N8

0

0

0

NR 0

0

0

Nx 0 0
0

ATi

0

0

Vyl

NI Vy

Vzl

N, vz

0

iVi

N*
(6.37)

式(6.37)における各節点の重み関数NAi ‑ 1

は,局所座標系における境界面の重心の座標vsff>VgiS>sを代入した式(6.33)である.

84

第7章 解析結果
本章では,今までの章で説明してきた手法による解析結果を説明する.最初に,描着剤の変形を記述するための材料モデルの検証を行う.材料モデルに関してはまず,粘弾性を記述する一般化Maxwellモデルの検証を行う.その後,ゴム弾性を記述できるよう修正したモデルの検証を行う.次に,移流計算テストの結果を説明し, 3次
元PLIC法による界面の捕捉精度を検証する.さらに,従来のEuler型固体達成手法では解析が困難であった,接触を伴う固体‑固体達成問題の解析結果を説明する.
最後に,固体境界面の‑軸引張問題の解析結果を用いて,従来のEuler型解法では取り扱うことのできなかった固体境界面に対する境界条件の付与が可能となったことを説明する.

7.1

一般化Maxwellモデルの検証

ここでは,本論文で取り扱う対象である粘着剤の変形を記述する材料モデルの検証を行う.まず,粘着剤の特徴の1つである粘弾性を記述できる材料モデルである
一般化Maxwellモデルについて,それを三次元に拡張し速度形で記述した構成方程式の妥当性を検証する.その検証には,応力緩和解析, ‑軸引張解析,自重変形解析という3つの解析を用いて行う.

7.1.1 応力緩和解析
最初の解析は応力緩和解析である.応力緩和は物体に一定ひずみを作用させている状態から時間経過と共に応力が減少してゆく現象であり,粘弾性材料における代表的な挙動である. Fig.7.1は解析モデルであり, Fig.7.2は境界条件を示している.
解析対象は長さが1cmの立方体であり,メッシュ分割数は, 3:方向, y方向, Z方向

85

にそれぞれ10である.扱う粘弾性体の密度は1000 kg/m3で微圧縮性を仮定するた
めポアソン比は0.49としている.

1.0 [cm]

1 .

O [cm]

1.0 [cm]

Fig. 7.1: Computational model for stress relaxation

Fig. 7.2: Boundary conditions for stress relaxation

本例題計算では一般化Maxwellモデルを7要素で表現する. Table 7.1は実験より得られた緩和時間と弾性係数である.初期状態において物体のZ方向に4.44×105 Paの垂直応力を与えるが,一般化Maxwellモデルの各要素について初期応力は,
El
(7.1)

0‑2 ‑4.44×105×

差Eiのように弾性係数の大きさに応じて与える.
Fig.7.3は厳密解と解析結果を比較したものである.なお解析解はFig.7.1における上端部の平均応力である.また厳密解は1次元の一般化Maxwellモデルの応力緩和を表現する以下の理論式[咽により求めている.

6‑差E*exp一芸)
86

(7.2)

Table7.1:

Relaxation

times

aノnd

Relaxation times回

elastic

moduliviscoelasticity

Elastic modul月Pa]

1.00 × lo一l0.1875 × 106

1.00× 100

0.1006 × 106

1.00 × 101

0.0427 × 106

1.00× 102

0.0326 × 106

1.00 × 103

0.0171 × 106

1.00 × 104

0.0100 × 106

1.00 × 1010

0.0530 × 106

50

of

100

150

Time [s]Fig. 7.3: Comparison between exact solution and computational result of stressrelaxation

Fig.7.3より応力が緩和することが再現されており,かつ厳密解と解析解が精度良く一致していることが確認できる.よって本計算において,本研究における構成方程式の応力緩和特性についてその正当性が確認できた.

87

7.1.2‑軸引張解析

次の解析は‑軸引張解析である. Fig.7.4は解析モデルであり, Fig.7.5はメッシュ
分割と境界条件を示している.物体は長さが1cmの立方体であり,オイラー型解析を行うために解析領域は物体よりも広く取っている.メッシュ分割数は, z方向, y
方向, Z方向にそれぞれ10である.粘弾性体の材料パラメータはすべて前述の応力緩和解析におけるものと同一である.

l.Ofcm]//1.0[cm]1.0 [cm

＼
2 .0 [cm ] ＼1.0 [c

X

2.
m ]

Fig. 7.4: Computational model for uniaxial tension

50.0 [mm/min]

50.0 [mm/mm]

蝣

+
賀蝣
H
M I
l
" M
l H
M I
l 至
=
…
i
…
=
】
鵜I
mm ,
E t i
註
≡
l
望
巨
盤出 轍 巨
I
I
"
表■
表■
表■
■■
一■
表■
■一
表■
‑
】
Ⅰ■
一題
■喜■
■
表■
一一
‑
旨
は確
固
野鶴
封宅

&

&

∴・・

∴

Fig. 7.5: Boundary conditions for uniaxial tension

有限要素法においては,節点にしか速度拘束を与えることができない.したがって物体が固定された計算メッシュを移動するオイラー型解法では,その物体端部に速度を作用させられない.よってここでの一章由引張解析においては, Fig.7.5に示す

ss

ように解析領域の上端部に外向きの速度を作用させて,メッシュ数は不変のまま等間隔に分割しなおすALE(Arbitrary Lagrangian‑Eulerian 型解法間を採用する.
Fig.7.6は本解析により得られた応力ーひずみ曲線と厳密解の比較である.解析解の応力は速度を与えた面の平均応力であり,厳密解は次式の1次元一般化Maxwell
モデルの構成方程式国から求めている.
ng・云T霊‑ ∑Elrl芸
:=1

(7.3)

i=l

Fig.7.6の両者はよく一致しており,構成方程式の引張特性とオイラー型解法が有効に動作していることが確認できる.

0%

1 0%

20%

30%

40%

50%

Nominal strainf‑1

Fig. 7.6: Comparison between exact solution and computational result of uniaxial
tension

7.1.3 自重変形解析
最後の解析はシリコーンパテの自重変形解析である.シリコーンパテの自重による変形は,長時間での流動でありその材料特性として粘弾性が支配的である.よって本論文における一般化Maxwellモデルによる構成方程式でも,十分な精度でシリ
コーンパテの自重変形解析が実施できるものと考えられる.
Fig.7.7は解析対象のシリコーンパテの大きざであり,解析のメッシュ分割と境界条件はFig.7.8に示す.,材料の密度は1095 kg/m3でポアソン比は微圧縮性とするため0.49としている.一般化Maxwellモデルの要素数はここでは6として, Table7.2
89

は実験より得られた緩和時間と弾性係数である.なお,本解析は自重変形解析であるので重力加速度は9.81 m/s2を作用させている.
五 O [cm ]l■" ■ ■ ら 一〇 〇〇 l 1.6[cm]1.5[cm ]2 .5
l.O rcm ]2.5[cm]Fig. 7.7: Computational model for silicone pate with deadweight:: : : :蝣蝣蝣蝣蝣
M ill 一
蝣
蝣蝣蝣蝣I
蝣
蝣蝣H I
蝣
HI
蝣
蝣蝣蝣蝣I
ほ: : : :IS !SB! i

&
l E I

!!!蝣!
!

&
M l
!蝣!

蝣
蝣蝣蝣蝣蝣蝣蝣蝣蝣I
i
i

I

&

&

Fig. 7.8: Boundary conditions for silicone pate with deadweightFig.7.9は実験結果と解析結果の比較である.それぞれの時刻における解析結果右側のコンターはミゼス応力[Pa]を示している.実験結果と解析結果とを比較して,定性的にではあるがその整合性が確認できる.なおこの解析は大変形を伴うため通常のラグランジェ型解法では計算メッシュが歪みすぎて解析が困難であるが,本論文におけるオイラー型解法を用いることで安定な解析が可能となる.

90

Table7.2: Relaxation times and elastic moduli of silicone pate

Relaxation times [s]elastic moduli [Pa]1.00 × 10‑2

2.15 × 108

1.00 × lo一l

3.22 × 105

1.00 × ‑10‑

5.34 × 103

1.00 × 101

3.46 × 102

1.00 × 102

9.41 × 100

1.00 × 103

7.12× 101

Fig. 7.9: Comparison between experiments and computational results of silicone
pate with deadweight

7.2

修正した一般化Maxwellモデルの検証

続いて,粘着剤の物性が示すもう1つの特徴であるゴム弾性を記述できる材料モ
デルの検証を行う.材料モデルは4.3節で示した一般化Maxwellモデルの修正モデ
ルである.

91

7.2.1 材料パラメータの取得
ここで,本論文で新たに提案した修正一般化Maxwellモデルの構成方程式(4.55)
を再掲する.
cr'

差2GiD+B{芸(T*‑ 2(1 +vl)Tl

抵差〈

3(1 ‑ 21/*)∴ ‑ KvtiD∫ (7.4)

この構成方程式を実際に数値シミュレーションで用いるためには,構成方程式内の材料パラメータを決定する必要がある.ここで扱っている構成方程式内の材料パラ
メータは, AiとBiそして7‑iである.これらの材料パラメータの同定は, ‑軸引張試験における実際の測定値と計算解とのカーブフィッティングにて行う.
ここで決意すべき材料パラメータを含む構成方程式(7.4)は3次元状態のものであるが,材料パラメータの同定で‑軸引張試験と整合させるために式(7.4)を1次元状態に縮約する必要がある.そこで式(7.4)を1次元状態に縮約すると次式のようになる.

雷‑差(E霊+B霊Ji一票(7.5)

式(7.5)は速度形であるが,陽的に時間積分することで任意のひずみ量における応力を計算することができる.
本論文における構成方程式は,修正した一般化Maxwellモデルの模式図(Fig.4.2)
からも分かるように,それぞれの要素で異なる緩和時間Tiにおいて対応する材料パ
ラメータAiとBiが必要になる. 4.3.1節でも述べたように,ここで扱う材料モデルの要素数は5要素であるからそれに対応する5種類の緩和時間Tiは1.00, 1.00× 101,
1.00×102, 1.00×10d, 1.00×106 を設定する.そしてこれら5種類の緩和時間に対応する材料パラメータAiとBiを同定する必要がある.緩和時間が5種類である
ので, AiとBiもそれぞれ5種類となり同定すべきパラメータは計10個となる.
これら計10個の材料パラメータを同定するための測定値としては, ‑軸引張試験
で10, 50, 300 mm/minの引張速度を与えて得た3種類の応カーひずみ関係を用いる.これら3種類の測定値と最も整合性が高くなるような10個の材料パラメータを,式(7.5)を用いて非線形最小二乗法により求める. Table.7.3に,最終的に得られた材料パラメータを示す.
92

Table7.3: Material parameters

Relaxation times [s]7.2.2

Al [Pa]B%

1.00 × 100

1.28 × 105

1.17

1.00× 101

8.93×101

4.40

1.00 × 102

3.71×102

3.67

1.00× 103

1.30×105

0.67

1.00 × 106

7.13×1013.85‑軸引張試験

ここでは実験結果との定量的整合性を検討するために,粘着剤の‑事由引張解析を扱う. Fig.7.10はその解析モデルである.解析対象は1辺1cmの立方体の粘着剤であり, Fig.4.2で示したような材料モデルの要素数は5として,材料パラメータは前節で同定したTable.7.3の値を用いる.また密度は1000kg/m3で,微圧縮性を仮定するために材料モデルのすべての要素でのポアソン比は0.49と設定する.本解析の目的は粘着剤の構成方程式の検証であるので,ここでは簡単な1格子での有限要素解析を行う.なお有限要素法における計算格子制御は,固体解析で一般的な計算格
子が物体の変形に追従するラグランジェ型解法を採用する.そしてFig.7.11は境界条件であり,物体変形の拘束を受けない均一変形となるような境界条件として,物体上端部のZ方向に強制速度を与える.

Fig. 7.10: Computational model for uniaxial tension

93Fig. 7.ll: Boundary conditions for uniaxial tension

解析の引張速度は5, 50, 500 mm/minとして,前節の材料パラメータ同定の際に用いた速度とはあえて3つのうち2つを違うものにする. Fig.7.12は引張速度5mm/minにおける実験結果と解析結果の比較であるがよく一致していることが確認できる.そしてFig.7.13とFig.7.14は引張速度がそれぞれ50 mm/minと500 mm/min
における実験結果と解析結果の比較であるが,ここでも両方ともで実験結果と解析
結果との良好な定量的整合性が確認できる.

200%

4 00%

600%

Nominal strata

Fig. 7. 12: Comparison between experiment and computational result in tensile speed
5mm/m上n

94

‑J
cd
'・蝣

邑
く/)

敬
Ql

.臼
こね

pa
s

・ ▼・・4

ME!
0.3

,*・

CI

0.2

2こ

0.1
0

200%

400%

600%

800%

Nominal strain

Fig. 7.13: Comparison between experiment and computational result in tensile speed
50 mm/min

0.8
I

I

氏.
'・ォ

0.7

‑E xperim ent
蝣C om putation

邑 0.6
【乃
亡′〉
01

上!
Zn

n

百

・I‑I

0.5
0.4
0.3

s
C)
Kt

200%

400%

600%

800%

Nominal strain

Fig. 7. 14: Comparison between experiment and computational result in tensile speed
500 mm/min

Fig.7.12, Fig.7.13, Fig.7.14の引張速度の異なるすべての解析において,ゴム弾性特有の応力の立ち上がりが再現できている.さらにFig.7.12とFig.7.13さらにFig.7.14とを比較すると,粘着剤のゴム弾性以外のもうひとつの特徴である粘弾性95

に関連する速度依存性も表現できていることが分かる.このことから本研究で提唱する新しい材料モデルは粘着剤の変形を記述する材料モデルとして妥当であることが判る.

7.3 移流計算テスト
本節では, Euler型解法において物質界面を捕捉する手法として検討した, 3次元PLIC法による界面捕捉の精度について, 2次精度風上差分法を用いたVOF法と比較しながら検証する.検証には, 2次元移流計算テストの結果と3次元移流計算テストの結果の合計3種類の解析結果を用いて示す.

7.3.1

2D Zalesak's disk rotation

Fig. 7.15: Geometry of Zalesak's disk prob‑
Fig. 7.16: Velocity丘eldlem

ここでは, 2D Zalesak's disk rotation問題[88]による移流計算テストの結果を示す.解析モデルはFig.7.15に示すように,半径0.15の円に幅0.05の溝がある形状である.速度場は,次式を与える.

‑‑My‑0.5)(7.6)
Vy‑+2tt(￡

96

0.5)

式(7.6)より1秒で図形が1回転することがわかる.速度場を図化したものをFig.7.16に示す.また,クーラン数は0.25,格子分割数は100× 100で計算を行った.
PLIC法による結果と2次精度風上差分法を用いたVOF法による計算結果をFig.7.20に示す. Fig.7.20より, VOF法では1回転後,境界面の数値拡散により図形の溝の形状を維持できなくなっているのに対して, PLIC法では急峻な境界面も捕捉できていることがわかる.
さらに,格子分割数642, 1282, 2562, 5122で計算を行い, Ll誤差により精度を定量的に比較した.またLl誤差の収束性を検討した. Ll誤差は,次式で定義している.

Ll‑志云腕‑￠￡71=1

ここで, Nは格子分割数>(7.7)^は初期のVOF値, 4,月ま最終的なVOF値(この問題

ではま‑1.00におけるVOF値)である. Fig.7.23に計算結果を示す.グラフの横軸
は格子幅の対数,縦軸はLl誤差の対数である. Fig.7.23より, VOF法では0.75次精度であるのに対して, PLIC法では1.12次精度で収束していることが確認できる.
また,グラフ中の白丸はVOF法におけるLl誤差の大きさに対するPLIC法におけるLl誤差の大きさの比である. N‑642では0.13, TV‑1282では0.07, JV‑2562
ではO.07, N‑h122では0.16となっており, PLIC法では誤差が平均して約92%減
少している.

7.3.2

2D Shearing now

ここでは, 2D Shearing flow問題[78]による移流計算テストの結果を示す.解析
モデルはFig.7.17に示すように,半径5/7T,中心座標>/2,0.2+tt/5)の円である.
速度場として,まず次式を与えて円をせん断変形させる. Fig.7.18は以下の速度場を図化したものである.
vz ‑ +sinxcosy(7.8)
Vy ‑ ‑cosxsin2/

さらに, 1000ステップ後,式(7.8)の逆向きの速度場,すなわちvx ‑ ‑sinxcosy(7.9)i) ‑+cosxsmy
97

Fig. 7.17: Geometry of shearing 凸ow problem Fig. 7.18: Velocity丘eldを与え,変形した図形を初期形状に戻す.クーラン数は0.25,格子分割数は100 ×100である.
PLIC法による結果と2次精度風上差分法を用いたVOF法による計算結果をFig.7.21に示す. Fig.7.21より, VOF法では界面が細くなるにつれて数値拡散が大きくなり,
最終形状は初期形状と大きく異なっている.一方, PLIC法では1000ステップ時には界面の分裂が見られるものの,細い界面形状も捕捉されており,最終形状がほぼ円形に回復していることが確認できる.
さらに,格子分割数642, 1282, 2562, 512^で計算を行い, Ll誤差により精度を定量的に比較した.またLl誤差の収束性を検討した. Ll誤差は,式(7.7)で定義している. Fig.7.24に計算結果を示す.グラフの横軸は格子幅の対数,縦軸はLl誤差の対数である. Fig.7.24より, VOF法では0.74次精度であるのに対して, PLIC法
ではO.95次精度で収束していることが確認できる.また,グラフ中の白丸はVOF
法におけるLl誤差の大きさに対するPLIC法におけるLl誤差の大きさの比である.
7V‑642では0.22, N‑ 1282では0.17, N‑2562では0.14, N‑5122では0.15となっており, PLIC法では誤差が平均して約83%減少している.

7.3.3

3D Stan ford bunny rotation

ここでは, 3次元複雑形状に対するPuC法の界面捕捉精度を検討する. 3次元複雑形状として, Fig.7.19に示すStanford bunny[91]を取り上げる.初期の物体形状
98Fig. 7.19: The stan ford bunnyについては, PLY形式のデータファイルをV‑Xgen[92]を用いてボクセルデータに変換したものを利用している.
解析領域の大きさは0.5×0.5 ×0.5であり,次式の速度場を与える.
V。 ‑ ‑2ir{y ‑yc)
Vy ‑ +2tt(∬‑ zc)(7.10)vz‑0

式(7.10)において x,とycはそれぞれ解析領域の中心点の￡座標とy座標で, zc‑‑0.01684, γ。ニー0.05である.またクーラン数は0.25,格子分割数は2563である.
PLIC法による結果と2次精度風上差分法を用いたVOF法による計算結果をFig.7.22に示す Fig.7.22はVOF値0.5の等数値面を措いたものである. Fig.7.22より, 1回
転後, VOF法ではStanford bunnyの耳の部分のように長細い界面を捕捉できていない.一方, PLIC法では初期形状をほぼ保っており,複雑形状に対する捕捉精度の高さを確認できる.
さらに,格子分割数323, 643, 1283, 2563で計算を行い, Ll誤差により精度を定量的に比較した.またLl誤差の収束性を検討した. Ll誤差は,式(7.7)で定義している. Fig.7.25に計算結果を示す.グラフの横軸は格子幅の対数,縦軸はLl誤差の対数である. Fig.7.25より, VOF法では0.56次精度で収束しているのに対して, PLIC
法では0.86次精度で収束していることがわかる.また,グラフ中の白丸はVOF法におけるLl誤差の大きさに対するPLIC法におけるLl誤差の大きさの比である.
99

7V‑323では0.41, TV‑643では0.33, N‑1283では0.26, iV‑2563では0.22となっており, PLIC法では誤差が平均して約70%減少している.

100

t‑0.00 (PLIC)

t‑0.00 (VOF)

t‑0.25 (PLIC)

t‑0.25 (VOF)

t‑0.50 (PI」IC)

t‑0.50 (VOF)

t‑0.75 (PLIC)

t‑0.75 (VOF)

t‑1.00 (PLIC)

t‑1.00 (VOF)

Fig. 7.20: Result of Zalesak's disk rotation

0 step (PLIC)

0 step (VOF)

500 step (PLIC)

500 step (VOF)

1000 step (PLIC)

1000 step (VOF)

1500 step (PLIC)

1500 step (VOF)

2000 step (PIJIC)

2000 step (VOF)

Fig. 7.21: Result of 2D shearing 丘ow

l

∴

1

i
I

i

I
I

I

i

L̲
t‑0.00 (PLIC)

t‑0.00 (VOF)a」一‑̲̲̲̲‑t‑0.25 (PLIC)ー̲

t‑0.25 (VOF)ヽ
t‑0.50 (VOF)P
t‑0.75 (VOF)

t‑1.00 (PLIC)

t‑1.00 (VOF)

Fig. 7.22: Result of 3D Stan ford bunny rotation

l

0.16

0.14

0.12

0.10

・‑7

(a

三・、

0.08

もO
CI

且̲9

0.06

0.04

0.02

0.00

‑10

‑9

‑8

‑7

‑6

‑5

x‑log2(h) (h: mesh size)

Fig. 7.23: L¥ error of Zalesak s disk rotation

0.32

0.28

0.24

0.20
o
:蝣サ
O

0.16

J
rLJ
bL)
CI

‑i

0.12

0.08

0.04

0.00
‑10

‑9

‑8

‑7

‑6

x‑log2(h) (h: mesh size)

Fig. 7.24: Lx error of 2D shearing flow

104

‑5

o
ヒ
O

J
N
CJj
〇

五‑ 0
‑ll

̲12

‑13

‑9

‑8

‑7

‑6

‑5

‑4

x‑log2(h) (h: mesh size)

Fig. 7.25: L¥ error of 3D stan ford bunny rotation

7.3.4 計算時間
次に, PLIC法とVOF法の計算時間を比較する. Fig.7.26‑Fig.7.28に,それぞれの移流計算テストに要した計算時間(CPU time)を示す.またFig.7.26‑Fig.7.28のグラフにおける白丸は, VOF法の計算時間を1とした場合のPLIC法の計算時間である. Fig.7.26‑Fig.7.28より,全ての移流計算テストにおいてPLIC法の計算時間の方が少ないことがわかる.具体的には, VOF法と比べPLIC法の計算時間は,2次元問題では平均して約20%程度少ない.
3次元PLIC法のアルゴリズムの実装に関しては,計算時間が非常に長くなることを予想していたが,上述したように実際には計算時間が減少するという結果が得られた.これには次の2つの理由が考えられる.第一に,境界面の法線ベクトルの計算について, 5.2.3節で述べたように,アルゴリズムが単純で計算コストの軽いYoungsの方法を用いたことである.移流計算の際, Youngsの方法を用いればステンシルは3格子でよいが, 2次精度風上差分法を用いたVOF法では5格子分のステンシルが必要になる.第二に,境界面の位置を決定する際に計算する3次方程式の解法としニフて,反復計算の必要がない方法を用いたことである.この手法では場合分けによって3次方程式の厳密解を直接的に計算するため,処理が非常に単純である.

105

53
∽
l

8

l

^^
qJ

.員
」6
P‑i
U
iさtちi
rヽ1

g)4

‑10

‑9

ー7

Iog2(h) (h: mesh size)

Fig. 7.26: CPU time of 2D Zalesak's disk rotation

‑10

‑9

ー7

Iog2(h) (h: mesh size)

Fig. 7.27: CPU time of 2D shearing 免ow

106

o

且10
LIT
4)

.g

二、
0,
U
DO
0

6

‑6

Iog2(h) (h: mesh size)

Fig. 7.28: CPU time of 3D Stan ford bunny rotation

1077.4 固体一固体達成問題
voF法を用いた界面捕捉法では国体の接触を伴う問題では解析結果は不合理な挙動となっていた.本節では,固体の接触を伴う問題に関して界面捕捉法としてPLIC法を用いた場合の解析結果について説明する.

7.4.1 2個の超弾性体の自由落下問題
ここでは,空気中で2個のneo‑Hooke体が自由落下する問題を解析する

neo‑

Hooke体の構成式に関しては後述の付録Aにて示す.解析モデルと物性値はそれぞれFig.7.29とTable.7.4に示す通りである.固体は弾性ゴムの1/100倍程度のせん
断弾性係数囲,流体は25‑Cの空気の物性値図であり,固体と流体で密度差は約1000倍ある.格子分割数は128×192である.また,この間題では計算の安定化のために1次風上差分法により速度の移流計算を行っている.
解析結果をFig.7.30とFig.7.31に示す.界面捕捉法にVOF法を用いた場合の解析結果も合わせて示している. Fig.7.30のt ‑ 0.46[s]で2個の固体が接触する前までは, PLIC法による結果とVOF法による結果は同じ挙動を示している.しかしt ‑ 0.76[s]以降においては, VOF法の場合,国体の接触により界面の数値拡散が激しくなり, 2個の国体は接着したような挙動を示し t‑2.00[s]においても2個の固体は接触したままの状態になっている.このように,接着したような挙動を示
すのは, VOF値の数値拡散により,基礎方程式が広範囲にわたって体積平均化されているためだと考えられる.一方, PLIC法の場合は境界面の数値拡散は激減して
おり,接着したような不合理な挙動は解消されている. t ‑ 0.46[s]で接触した後,
t‑0.76[s]では固体同士の反発力により跳ね上がる様子が再現されている.同様に,
t‑0.98[s]で接触した後

t‑1.20[s]でも跳ね上がっている.そして, i‑1.48sからt ‑ 2.00[s]にかけて固体が接触した後,底面で運動が平衡状態に達している.

108

Fig. 7.29: Computational model of two free falling solids

Table

7.4: Material properties

Solidl (neo‑Hookean model)

Solid density ββ

10叫kg/m3]

Shear modulus G 0.50 × 104[Pa]Solid2(neo‑Hookean model)

Solid density β3

10叫kg/m3]

Shear modulus G 1.00 × 104[Pa]Fluid(Newtonian臥Iid)
Fluid density pf 1.184[kg/m3
Viscosity fi 1.82 × 10"5[Pa‑s]

gravitational acceleration

109

9.叫m/S2]

t‑0.00 s] (PLIC)

t‑0.00回(VOF)

t‑0.32 s (PLIC)

t‑0.32 s (VOF)

t‑0.46 s (PLIC)

t‑0.46[s] (VOF)

t‑0.76 s (PLIC)

t‑0.76 s] (VOF)

Fig. 7.30: 酔ee fall of two neo‑Hookean solids (0 ≦ i ≦ 0.76)

t‑0.98回(PLIC)

t‑0.98 s (VOF)

t‑1.20[s] (PL工C)

t‑1.20 s (VOF)

t‑1.48s (PLIC)

t‑1.48 s (VOF)

t‑2.00 s (PLIC)

t‑2.00 s (VOF)

Fig. 7.31: nee fall of two neo‑Hookean solids (0.98 ≦ i ≦ 2.00)

7.5 固体境界面の一幸由引張問題
本節では,力学的境界条件及び幾何学的境界条件をPLIC法で定義した国体境罪面に付与することにより,国体の‑軸引張問題を解析する.

7.5.1 力学的境界条件による一軸引張
ここでは,力学的境界条件による‑軸引張問題の解析結果を説明する. Fig.7.32に示すように,固体境界面に分布荷重を作用させる.固体以外の領域は流体で満たされており, Ⅹ軸の正の方向に半無限に解析領域が連続しているものとし,奥行き方向の長さは6.25 × 10‑2である.固体と流体の物性値はTable.7.5の値を与える.
本解析の目的は構築した境界条件付与手法の妥当性を検証することであるため,物性値としては無次元化した数値を用いた.分布荷重は,固体の振動を避けるために,
Fig.7.33に示すように線形に増加させている.解析領域端部の境界条件はローラ支点とし,格子分割数は64×64である.
Fig.7.34は固体の変形の様子を示しており,力学的境界条件によって固体が上方に引っ張られている様子が確認される.
また, Fig.7.35に示すグラフは,数値解と理論解を比較したものである.縦軸が偏差応力のzz成分,横軸が固体の長さである. neo‑Hooke体の‑軸引張の理論解は,次式により与えられる囲.
0‑乞‑G(f‑1)(7.ll)

ここで, lは固体の長さである. Fig.7.35より,数値解と理論解は精度良く一致していることが確認できる.よって,従来のEuler型解法では困難であった力学的境界条件の付与が可能になったといえる.

112

load per unit

F lu id

D istri b u te d lo a d
2 ‑0

1

1

1

⊥

ユ

⊥

1

1

…1 .0
z
.lLL

ー
lLL

ー
lLL

lユ⊥

2 .0
X

Fig. 7.32: Model of uniaxial tensile Fig. 7.33: Time history of distributed load

Table

7.5: Material properties

Solid(neo‑Hookean model)
Soliddensity/フβ

1.0Shear modulus G 1.0
Fluid(Newtonian auid)
Fluid density p/
viscosity /j,

113

1.0
1.0 × 10‑2

t=0.0Fig. 7.34: Deformation of neo‑Hookean solid

8. OE‑02

7.OE̲02

6.OE‑02

5.OE‑02
GO
GO
O
上]

ロケ1
く)

'己
CI
cd
>・
O
'13

4.OE‑02

3.OE‑02

2. OE‑02

1.OE‑02

O.OE+00

1.00

1.05

1.10

1.15

1.20

1.25

1.30

lengthFig. 7.35: Comparison of deviatric stress zz variation with exact solution

114

7.5.2 幾何学的境界条件による‑軸引張
ここでは,幾何学的境界条件による‑軸引張問題の解析結果を説明する. Fig.7.36に示すように,固体境界面を速度vz‑3.0×10 3を与えて引っ張る.固体以外の領域は流体で満たされており, Ⅹ軸の正の方向に半無限に解析領域が連続しているものとし,奥行き方向の長さは6.25 × 10‑2である.固体と流体の物性値はTable.7.5の値を与える.解析領域端部の境界条件はローラ支点とし,格子分割数は32×32である.

Fig. 7.36: Model of uniaxial tensile with geometrical boundary condition

Fig.7.37は固体の変形を示しており,幾何学的境界条件によって固体が上方に引っ張られている様子が確認される.また, Fig.7.38に示すグラフは,数値解と理論解を比較したものである.縦軸が偏差応力のzz成分,横軸が固体の長さである. neo‑Hooke
体の‑軸引張の理論解は,前述したように式(7.ll)により与えられる. Fig.7.38より,
数値解と理論解は精度良く一致していることが確認できる.よって,従来のEuler
型解法では困難であった幾何学境界条件の付与が可能になったといえる.

115

蝣*

* 'm

m'

m蝣

K

,

^

'

蝣
蝣
蝣
‑蝣
^ ‑"4 ':

t=0.0

Fig. 7.37: Deformation of neo‑Hookean solid

め
帆
V
上]

3 0.04‑⊂ocd

̀己0.03
'tコ

1.00

1.05

1.10

1.15

1.20

1.25

1.30

lengthFig. 7.38: Comparison of deviatric stress zz variation with exact solution

116

第8章 結言
本章では,今までに述べてきた粘着剤の大変形解析に関する研究をまとめる.
本研究の目的は,大変形を伴う粘着剤の変形解析を可能とすることである.粘着剤の変形解析を行うには以下の課題があった.
1.粘着剤の変形は大変形を伴うため,固体解析にて通常用いられるLagrange型
解法では計算格子の破綻により,解析が困難である.
2.粘着剤はゴム弾性及び粘弾性という性質を有しており,その変形を記述することが難しい.
3. Euler型解法の使用においては,長時間積分においては数値拡散が生じ,また物質界面に対して境界条件を付与しにくい.

大変形解析に対してはEuler型解法を適応することで解析を可能とした.粘着剤の変形の記述に関しては新しい材料モデルを提唱した.物質界面への境界条件付与に関しては, PLIC法を用いた物質界面の高精度化を行った.
以下に各章で得られた要点をまとめて,本論文の結言とする.

第1章では,本研究の目的と背景について述べた.粘着テープは様々な産業分野で用いられている製品であるが,その設計に数値計算は殆ど用いられていない.その理由は粘着テープの構成要素である粘着剤 が容易に大変形すること,及びその変形挙動が非常に複雑なためであると述べた.今後,粘着テープに対してさらに高機能化が要求されるなか,粘着剤の設計に対しても数値解析による設計支援が必要不可欠となる.よって,本研究の目的は非線形・大変形を示す粘着剤の変形を記述できる数値解析手法を構築することである.

117

第2章では,粘着剤の力学特性について述べた.粘着剤は高分子で構成されているため,その弾性は高分子の絡まり合いに起因するエントロピー変化により生じる
ゴム弾性であること,さらに,粘性と弾性を併せもつ粘弾性体であり応力は時間依存性を有していることを述べた.

第3章では,粘着剤の大変形解析を可能とするために用いたEuler型有限要素汰
について述べた.まず,有限要素法の基礎方程式である保存則について説明した後,Euler記述の方程式を解くための手法であるOperator Split法を説明した. Operator
Split法の非移流ステップにおけるLagrange型解法に関して,本研究で用いている
動的陽解法を説明した.また,移流ステップで用いる界面捕捉法についてVOF法を説明した.

第4章では,粘着剤の複雑な変形を記述する構成方程式について述べた.まず,粘着剤の有する粘弾性を記述するために,材料モデルとして一般化Maxwellモデルを用いた.そして一般化Maxwellモデルを三次元に拡張し,その構成方程式について
速度形の定式化を行った.この際,客観応力速度としてJaumann速度を用いた.衣いで,粘着剤が有するもう1つの特性であるゴム弾性を記述するために,先に定式化した一般化Maxellモデルのスプリング部を定数からひずみの関数とするよう修正した新しい材料モデルを提案した.新しい材料モデルにおいても三次元に拡張した速度形の構成方程式を導いた.この新しい材料モデルを用いることで,粘着剤が示す非線形・大変形,かつ速度依存性を有する複雑な変形挙動を記述することができる.

第5章では, Euler型有限要素法を用いることで生じた課題である長時間積分における数値拡散に関して,その解決手法として採用した高精度な界面捕捉法としてのPLIC法を説明した.本研究では三次元に拡張したPLIC法を用いた.その際,界面法線ベクトルの計算には非移流ステップで求めるVOF関数を用いた.また,境界セルの距離定数を求める際に生じる三次関数の計算では,収束計算を用いず直接計算する手法を採用した.このようにして構成した界面捕捉法によって, Euler型解法においても高精度に物質界面を定義することが可能となる.

118

第6章では第5章にて述べたPLIC法を用いて陽に定義された物質界面に,境界条件を付与する方法について述べた.力学的境界条件は仮想仕事式を導く際に自動的に取り込まれるが,幾何学的境界条件は取り込まれない.本研究では幾何学的境界条件については,ペナルティ法を用いて仮想仕事式に取り入れる手法を採用した.
こうすることで,従来のEuler型有限要素法では困難であった,物質境界面への境界条件付与が可能となる.

第7章では前章までに述べた手法の妥当性を検証するために行った解析について述べた.粘弾性を記述するための材料モデルとして用いた一般化Maxwellモデルの検証においては,応力緩和試験, ‑事由引張試験,自重変形解析による解析を行い,実験結果とよい一致をすることを示した.このことにより,一般化Maxwellモデルを三次元へ拡張したことの妥当性を検証できた.次いで,ゴム弾性を示すために一般化Maxwellモデルを修正した新しい材料モデルの妥当性を検証するために, ‑事由引張試験を用いた解析を行った.結果,ゴム弾性の指標となる‑車由伸張変形の応カーひずみ曲線において計算結果と実験値との間でよい一致を示した.解析は単一の引張速度だけでなく複数の引張速度でも行い,非線形大変形を示す応力‑ひずみ曲線の時間依存性についても解析結果と実験値でよい一致を示した.このことから,新しい材料モデルを用いることにより,粘弾性・ゴム弾性という2つの性質を有する粘着剤の複雑な変形を記述できるようになったことを示した.
また, PLIC法による界面捕捉精度の向上に関して,二次元問題と三次元問題について解析を行った.結果,従来の界面捕捉法であるVOF法に比べてPLIC法を用いた手法では極めて高精度に界面を保持していることを示した.さらに,解析時間についてもPLIC法を用いた方が早くなることを示し,本手法の有効性を示した.そして, PLIC法により定義された物質界面に対する境界条件の付与法を検証するために,力学的境界条件と幾何学的境界条件の各々について‑軸引張解析を行った.結栄,計算解と理論解との間でよい一致をすることが確認され,本手法の妥当性を示すことができた.

以上より,本研究の目的は達成されたと考える.今後は,粘着剤の実際の評価に適応することが課題となる.具体的な課題としては,以下の5点が考えられる.

119

1.粘着剤は微圧縮性を示すが,本研究における構成方程式の体積成分では線形圧
縮性弾性体の構成式を用いている.そのため,数値解析が不安定になる可能性がある.微圧縮性を表現するために圧力を精度良く記述できる手法の検討が必要である.
2.実際の粘着剤の評価を解析する場合,たとえばピール試験では,その際に生じる糸曳きを解析しなければならない.この糸曳きには剥離や破壊といった現象が生じる.これら現象を記述できる手法の検討が必要である.
3.実際の粘着剤は物性や変形に温度依存性がある.温度一時間換算則18]などを利用して温度依存性を記述する必要がある.
4.本手法では,速度を算出する際に陽解法を用いており,クーラン条件により時間増分が制限されるため,時間スケールの長い間題や固体の弾性係数が大きい問題を解析するには莫大な計算時間が必要となる.よって,時間増分を大きく設定できるNewmarkβ法などの陰解法を検討する必要がある.
5.現状の手法では,境界条件を与える位置をLagrange的に追跡しないため,境界条件を与える位置の制御が難しい.よって,境界条件を課す位置を任意に制御できる手法を導入する必要がある.

120

付録A

非圧縮性固体及び流体の構成
方程式

本章では, 7章で用いる非圧縮性超弾怪体の‑種である非圧縮性neo‑Hooke体の構成方程式について説明する[95]‑[97].また達成解析で用いる非圧縮性Newton流体の構成方程式についても説明する[57]‑64 95‑ 971.

A.1

超弾性体

A.1.1

非圧縮性neo‑Hooke体

超弾性体とは,次式のように,変形やひずみの成分によって微分されることにより共役な応力成分を生じる弾性ポテンシャル関数Wが存在する物質である.

S‑芸

(A.I)

ここで, Sは第2Piola‑Kirchhoff応力テンソル, EはGreen‑Lagrangeひずみテンソ
ルである. SとEは双方とも観測不変テンソルで, Wは客観性のあるスカラーであることより,超弾性体の構成方程式(A.I)は観測者によらず同一であり,物質客観性の原理を満たしている.本研究では,右Cauchy‑Green変形テンソルCひいては左Cauchy‑Green変形テンソルBを用いて構成方程式(A.I)の定式化を進めることにする.まず, E‑{C‑I)/2であることより,式(A.I)を次式のように変形する.

S‑ 2芸

(A.2)

一般に弾性ポテンシャル関数は,右Cauchy‑Green変形テンソルCの主不変量存,
He, Hidの関数として与えられることより,偏微分の連鎖律の公式を用いて,

5‑2芸雷+蒜等+蒜等)
121

(A.3)

と表すことができる. Cの主不変量のCに関する偏微分がそれぞれ
I

IcI‑C
IIIcC‑1となることを用いれば,式(A.3)は次式になる.

S‑2芸十品Ic)‑蒜C・蒜nice‑i (A.7)<t‑F‑S‑FT/J及び存‑IB, IIc‑Hb, HIc‑Mbが成り立つことより,

式(A.7)は

a‑J荒IIIbI+芸+品IB) B一芸B・B(A.8)と表せる.式(A.8)のように,超弾性体の構成方程式はCauchy応力Uと左Cauchy‑
Green変形テンソルBを用いて表せる.以上の定式化は圧縮性超弾性体に対するものであり,非圧縮性が仮定されれば,式(A.8)は次式のように修正される.すなわち,主不変量のうちHIb‑1となることから,弾性ポテンシャル関数WはIb, Hbのみの関数となること,及びJ‑1を考慮して次のようになる.

cr‑2芸+荒Ib)B‑2荒B B‑pI (A.9)
ここで, pは不定圧力であり,応力の体積成分は物質点の運動の履歴からは定める
ことはできず,境界条件により決定される.このように,非圧縮性物質では等積変形のみが可能であり,応力決定の原理は制約を受ける.
本研究で取り扱う非圧縮性neo‑Hooke体は,次式で定義される弾性ポテンシャル関数を持つ物質である.
W‑

ci(/c‑3)(A.10)‑

C!(/B‑3)(A.ll)

式(A.ll)において, clは実験により定められる定数であり,微小変形時にはHooke
則のせん断弾性係数Gとの間にcl ‑ G/2なる関係が成り立つ.非圧縮性neo‑Hooke
体の弾性ポテンシャル関数(A.ll)を非圧縮性超弾性体の構成方程式(A.9)に代入す
ることにより,非圧縮性neo‑Hooke体の構成方程式は次式になる.
a‑GB‑pI
122(A.12)

A.I.2

左Cauchy‑Green変形テンソルの評価

固体力学で一般に用いられるLagrange表示では,物質点の初期の位置ベクトルX
を用いて物質の変形を評価する.そのため,変形勾配テンソルFひいては左Cauchy‑
Green変形テンソルBを容易に求めることができる.これに対して,本手法ではEuler
表示により固体の変形を記述するため,変形勾配テンソルを直接的に求めることができない.よって,左Cauchy‑Green変形テンソルβの評価には工夫が必要となる.
そこで本研究では,次式に示すように左Cauchy‑Green変形テンソルβの定義式
B=F Flの両辺を物質時間微分することにより得られる恒等式を導入する.
DB(A.13)

D(F蝣F7

Di

Dt

DFT

芸?T+F‑‑㌃
L

F

Fl+F‑Fl蝣LI

L B+B LT(A.14)(A.15)(A.16)

ここで, Lは速度勾配テンソルである.さらに式(A.16)をEuler表示すれば,次式になる.

雷 v∇B‑LB+BLT(A.17)

本研究では,式(A.17)によりEuler表示の左Cauchy‑Green変形テンソルβを時間発展させる.

A.2

非圧縮性Newton流体

本研究では,水などの液体や音速に比べて流速が小さい気体を取り扱う.よって,流体については非圧縮性Newton流体の構成方程式を用いる. Newton流体とは,せん断応力がせん断変形に比例するという法則にしたがう流体であり,構成方程式は次式で与えられる.
a‑2/jD+<‑p+ I k‑‑/j, ItrD

(A.18)

ここで, FLは粘性係数, Dは変形速度テンソル, pは流体の圧力, Itは体積粘性率である.非圧縮性物質の場合, tiD=0が成り立つことにより,構成方程式は以下の123

ようになる.

g‑ 2/JD‑pZ

124(A.19)

謝辞
本論文は筆者が広島大学大学院工学研究科社会環境システム専攻博士後期課程に在籍中の研究成果をまとめたものです.本研究を進めるにあたっては,大変多くの方々の御世話になりました.ここに深く感謝の意を表します.
同専攻准教授岡津重信先生には,指導教員として研究活動全般にわたり大変に熱
心で温かいご指導を賜りました.論文の書き方のような基礎から,本研究の本質まで親切丁寧にご教授していただき深謝しています.同専攻教授北村充先生にはオハ
イオでの始めての国際会議にて大変お世話になりました.深く感謝の念を申し上げます.同専攻教授藤本由紀夫先生には粘着テープに関心を示していただき,興味深
いご指摘をしていただき感謝申し上げます.同専攻教授藤井堅先生には本研究に貴重なご意見をいただき感謝申し上げます.同専攻助教田中智行先生には様々な相談をしたときでも,丁寧にご教授していただき大変感謝しています.
株式会社日立製作所の中尾賢司氏,広島大学大学院在学の岩本崇氏,慶応義塾大
学大学院在学の横山睦女史には本研究に携わっていただきご尽力頂いたことに感謝
申し上げます.
日東電工株式会社の西口浩司氏には,彼が広島大学に在籍中から本研究に携わっていただき並々ならぬ努力と貢献をしていただいたことに,深く感謝の念を申し上げます.
日東電工株式会社の主幹研究員浅野猛氏,主幹研究員佐野建志氏には粘着剤の変形に関して未熟な私をご指導していただき大変感謝しています.
このような大変貴重な機会を与えて頂いた日東電工株式会社の執行役員表利彦博士,主幹研究員杉原保則氏,主幹研究員大西啓仁氏,主任研究員梅田道夫氏,主任研究員田中良和氏に感謝申し上げます.
最後に,本研究を続けることにあたり,何も文句を言わず直向に家庭を支えてくれた妻洋恵に感謝申し上げます.
以上の皆様方にここに重ねて厚く謝意を表し,謝辞といたします.
125

参考文献
ll]浦浜圭彬‥粘着テープ剥離のモルフォロジー,日東技報 29,No.2,pp.1‑ll, 1991
[2] K. Hino and H. Hashimoto: Morphological studies on the adhesion mechanismof pressure sensitive adhesives, Journal of Applied Polymer Science, 30, pp.3369‑
3376,1985[3] A. Zosel: Adhesive Failure and Deformation Behaviour of Polymers, The Journalof Adhesion, 30, pp.135‑149, 1989

[4]

A.

Zosel‥

Physical

Properties

andAdhesion

Performance

of

PSAs,

Adhesives

Age, 32, No.ll, pp.42‑47, 1989
[5] D. J. Benson and S. Okazawa: Contact in a multi‑material Eulerian finite element
formulation, Computer Methods in Applied Mechanics and Engineering, 193.
pp.4277‑4298, 2004.

囲岡揮重信,河口篤志,藤久保昌彦‥各種メッシュ制御における動的陽解法,応用力学論文集6, pp.151‑158, 2003.
[7] S. Okazawa, K. Kashiyama and Y. Kaneko: Eulerian formulation using stabi‑
lized finite ele甲ent method for large deformation solid dynamics, International
Journal for Numerical Methods in Engineering, 72, pp.1544‑1559, 2007.
[8] D. L. Youngs: Time‑dependent multi‑material 凸ow with large月Iuid distortion,

Numerical Methods for Fluid Dynamics, edited by K. W. Morton and M. J.
Baines, Academic Press, pp.273‑285, 1982.
[9] E. Aulisa, S. Manservisi, R. Scardovelli and S. Zaleski: Interface reconstruction
with least‑squares fit and split advection in threedimensional Cartesian geometry,
Journal of Computational Physics, 225, pp.2301‑2319, 2007.
126

[10] J. Hernandez, J. Lopez, P. Gomez, C. Zanzi and F. Faura: A new volume of
fluid method in three dimensions. Part I: Multidimensional advection method
withface一matchedflux

polyhedra,

International

Journalfor

Numerical

Methods

in Fluids, 58, pp.897‑921, 2008.
[11] P. Liovic, M. Rudman, J. L. Liow, D. Lakehal and D. Kothe: A 3D unsplit‑
advection volume tracking algorithm with planaritypreserving interface recon‑

struction, Computers &動ids, 35, pp.1011‑1032, 2006.
[12] J. Lopez, C. Zanzi, P. Gomez, F. Faura and J. Hernandez: A new volume offluid
method in three dimensions. Part II: Piecewise‑planar interface reconstruction

with cubic‑Bezier fit, International Journal for Numerical Methods in Fluids,
58, pp.923‑944, 2008.

[13] D. Lorstad and L. Fuchs: High‑order surface tension VOF‑model for 3D bubble
flows with high density ratio, Journal of Computational Physics, 200, pp.153‑
176,2004.

[14] G. H. Miller and P. Colella: A conservative three‑dimensional Eulerian method
for coupled solid一缶uid shock capturing, Journal of Computational Physics, 183,

pp.26‑82, 2002.

15] Y. Renardy and M. Renardy: PROST: a parabolic reconstruction of surface
tension for the volume‑of‑fluid method, Journal of Computational Physics, 183.
pp.400‑421, 2002.
[16] A. V. Pocius, "Adhesion and Adhesives Technology: An Introduction." Hanser‑
Gardner Publishers, Munich, Germany, 1997

囲村上謙古"やさしいレオロジー",産業図書, 1986
[18]村上謙吉:レオロジー基礎論,産業図茸1991[19] C. A. Dahlquist: Pressure Sensitive Adhesives, In R. L. Patrick (ed.), Treatiseon adhesion and adhesives, Marcel Dekker, New York, 2, pp. 219‑270, 1969127

[20] R. S. Rivlin: The e庁ective work of adhesion, Paint Technology, 9, pp.215‑216,
1944[21]畑敏雄,高分子化学, 4, 67, 1947
B. V. Deryagin, and N. A. Krotova, Doklady Akademii Nauk SSSR, 61, p849.
1948[23] J. J. Bikerman: Experiments on Peeling, Journal of Applied Polymer Science,
2, pp.216‑224, 1959
[24] J. L. Gardon: Peel adhesion. II. A theoretical Analysis, Journal of Applied
Polymer Science, 7, 643, 1963
25] D. H. Kealble: Theory and Analysis of Peel Adhesion: Bond Stresses andDistributions, Transactions of the Society of Rheology, 4, pp.45‑73, 1960

[26]斎藤隆則‥粘着の理論I,日本接着協会識21, pp.228‑239, 1985
[27]三木祥男:連続体力学による剥離理論の定式化,日本接着協会識24, pp.99‑105,
1988

囲A. N. Ge叫Adhesion and strength of viscoelastic solids. Is there a relationshipbetween adhesion snd bulk properties, Langmuir, 12, pp.4492‑4496, 1996

囲W. C. Dale, M. D. Paster, and J. K. Haynes: Mechanical properties of acrylicpressure sensitive adhesives and their relationships to industry standard testing,
Journal of Adhesion, 31, pp・1‑20, 1989

[30]浦浜圭彬: "タック"いろいろなアプローチ,接着 38,No.10,pp.452‑463, 1994[31] E. P. Chang: Viscoelastic windows of pressuresensitive adhesive, Journal of
Adhesion, 34, pp.189‑200, 1991
[32] W. O. David, T. Nipithakul, M. Nardin, J. Schultz, and K. Suchi va: In且uenceof Nonrubber Constituents on Tack of Natural Rubber. I. At Very Short Timesof Contact (Pendulum Test), Journal of Applied Polymer Science, 78, pp.1486‑
1494,2000
128

[33] T. Hata: Mechanisms of Adhesive Failure, Journal of Adhesion, 4, pp. 161‑170.
1972

H. Mizumachi: Theory of Tack of Pressure Sensitive Adhesive. I, Journal ofApplied Polymer Science, 30, pp. 2675‑2686, 1985

囲H. Mizumachi, and Y. Hatano: Theory of Tack of Pressure Sensitive Adhesive.
II., Journal of Applied Polymer Science, 37, pp. 3097‑3104, 1989
M. Mooney: A theory of large elastic deformation, Journal of Applied Physics,
11, pp.582‑592, 1940

岡R.W. Ogden: Large deformation isotropic elasticity, on the correlations oHhe
theory and experiment for compressible rubber‑like solids, The Royal Society of
London, A‑326, pp.565‑583, 1972[38] J.C. Simo: On a fully three‑dimensional丘nite‑strain viscoelastic damage model:

formulation and computational aspects, Computer Methods in Applied Median‑ics and Engineering, 60, pp.153‑173, 1987
[39] T. Sussman, and K. Bathe: A finete element formulation for nonlinear incom‑pressible elastic and inelastic analysis, Computers & Structures, 26, No.1/2,
pp.357‑409, 1987

[40] T. Yamada, and F.Kikuchi: A mixed丘nite element method for large deforma‑tion analysis of incompressible hyperelastic materials, Theoretical and Applied

Mechanics, 39, pp.61‑73, 1990
[41] M. Doi, and S. F. Edwards: The Theory of Polymer Dynamics, Clarendon
Press, Oxford, 1986[42] E. M. Arruda, and M. C. Boyce: A three‑dimensional constitutive model for thelarge stretch behavior of rubber elastic materials, Journal of Mechanical Physics
and Solids, 39, No.2, pp.389‑412, 1993129

[43] Y. Tomita, T. Adachi, and S. Tanaka: Modelling and application of constitutive
equation for glassy polymer based on nona:蝕Ie network theory, European Journalof Mechanics A/Solids, 16, No.5, pp.745‑755, 1997[44] Y. tomita, K. Azuma, and M. Naito: Strain‑rate‑ dependent deformation be‑havior of corbon‑black filled rubber under monotonic and cyclic straining, KeyEngineering Materials, 340‑341, pp.1017‑1024, 2007

囲A. S. Lodge: Elastic liquid, Academic Press (London), 1964
[46]三木祥乳浅野猛‥粘弾性体の大変形に関する有限要素法による定式化,日本接着協会誌, 22, pp.463‑472, 1986
[47]

B.

C.

plastic

Duncan‥FailureMaterials,

Criteria

Project

PAJl‥

and

their

Final

Application

Report,toNational

Visco‑elastic/Visco‑Ph

ysicalLaboratory,

Teddington, UK, Octobe, 1999.

[48]朱峰承興:粘着剥離挙動と粘着伸張変形との相関,日本接着学会識43, pp.481‑
487,2007

T. Yamaguchi, H. Morita, and M. Doi: Modeling on debonding dynamics ofpressure‑sensitive adhesives, European Physical Journal E, 20, pp.7‑17, 2006[50] http://octa.jp, OCTA Home page.
[51] H. Motrita, M. Yamada, T. Yamaguchi, and M. Doi: Molecular dynamics studyof the adhesion between end‑grafted polymer films, Polymer Journal, 37, No.10,
pp.782‑788, 2005.

囲藤井幸蔵‥流体力学の数値計算法,東京大学出版会, 1994.
B. van Leer: Towards the ultimate conservative difference scheme IV:a new ap‑
proach to numerical convection, Journal of Computational Physics, 23, pp.276‑
299,1977.

130

[54] B. van Leer: Towards the ultimate conservative difference scheme V:a second‑order sequel to godunov's method, Journal of Computational Physics, 32,
pp.101‑136, 1979.

55] C. W. Hirt and B. D. Nichols: Volume offluid (VOF) method for the dynamics
of free boundaries, Journal of Computational Physics, 39, pp.205‑26, 1981.

囲数値流体力学編集委員全編, ‥数値流体力学シリーズ4 移動境界流れ解析東京大学出版会, 1995
岡久田俊明:非線形有限要素法のためのテンソル解析の基礎,丸善, 1992.
[58] G. A. Holzapfel: Nonlinear solid mechanics / a continuum approach for engi‑

neering, John Wiley and Sons, 2000.
[59] J. Bonet and R. D. Wood: Nonlinear continuum mechanics for finite element
analysis (2nd Edition), Cambridge University Press, 2008.
[60] A. J. M. Spencer: Continuum mechanics, Dover publications, 1980.
[61] J. E. Marsden and T. J. R. Hughes: Mathematical foundations of elasticity,
Dover publications, 1983.

囲京谷孝史‥よくわかる連続体力学ノート,森北出版, 2008.
[63]棚橋隆彦:連続体の力学(1) ‑物質の変形と流動‑ (第3版) ,理工図書, 1995.
[64]石原繁:テンソル‑科学技術のために‑ (第12版) ,裳華最2008.
[65] A. J. Chorin: Flame advection and propagation algorithm, Journal of Compu‑
tational Physics, 35, pp.1‑31, 1980.

囲久田俊明,野口裕久:非線形有限要素法の基礎と応用,丸善, 1995.
[67]日本塑性加工学会編:非線形有限要素法一線形弾性解析から塑性加工解析まで‑,コロナ社, 1994
[68] A. Harten: On a class high resolution total‑variation‑stable finite difference
schemes, SIAM Journal on Numerical Analysis, 29, pp.329‑349, 1981.
131

[69]日本計算工学会流れの有限要素法研究委員全編:続・有限要素法による流れのシミュレーション,シュプリンガ‑・ジャパン, 2008.

囲山本三三三‥材料科学のための物体の変形学,誠文堂新光社, 1972
[71]富田幸雄:レオロジーー非線形流体の力学‑,コロナ, 1975
囲横道英雄:テンソルとレオロジー,技報堂, 1975
囲山田嘉昭:塑性・粘弾性,培風館1980[74] P. Chadwick :連続体力学一簡明な理論と例題‑,ブレイン図書出版, 1979[75] D. J. Benson: Computational methods in Lagrangian and Eulerian hydrocodes,

Computer Methods in Applied Mechanics and Engineering, 99, pp.235‑394, 1992.
[76] W. J. Rider and D. B. Kothe: Reconstructing volume tracking, Journal ofComputational Physics, 141, pp.112‑152, 1998.
[77] R. Scardovelli and S. Zaleski: Direct numerical simulation of丘ee‑surface andinterfacial flow, Annual Review of Fluid Mechanics, 31, pp.567‑603, 1999.
[78] M. Rudman: Volume tracking methods for inter facial flow calculations, Inter‑national Journal for Numerical Methods in Fluids, 24, pp.671‑691, 1997.
[79] R. Scardovelli and S. Zaleski: Analytical relations connecting linear interfaces
and volume fractions in rectangular grids, Journal of Computational Physics,
164, pp.228‑237, 2000.
[80] A. J. Chorin: Curvature and solidification, Journal of Computational Physics.
58, pp.472‑490,1985.
[81] T. J. Barth: Aspects of unstructured grids and finite‑volume solvers for Eulerand navier‑stokes equations, VKI/NASA/AGARD Special Course on Unstruc‑tured Grid Methods for Advection Dominated Flows AGARD Publication R‑787,
1995.

132

[82] B. Swartz: The second‑order sharpening of blurred smooth borders, Mathemat‑
ics of Computation, 52, pp.675‑714, 1989.

M. Metcalf, W. H. Press, S. A. Teukolsky, W. T. Vetterling and B. P. Flannery
Numerical recipes in Fortran 90, Cambridge University Press, 1996.

[84] D. J. Benson: Volume o用uid interface reconstruction methods for multi一material problems, Applied Mechanics Reviews, 55, pp.151‑165, 2002.
[85] M. Abramowitz and I. A. Stegun: Handbook of Mathematical Functions, Dover
Publications, 1965.

[86]長嶋利夫,新山健二,石原嘉一:構造的なメッシュを用いた応力解析における基本境界条件付与方法の検討,日本機械学会論文集(A編 70, pp.354‑362,
2004.

岡T. Belytschko, W. K. Liu and B. Moran: Nonlinear丘nite elements for continua
and structures, John Wiley and Sons, 2000.
S. T. Zalesak: Fully Multidimensional Flux‑Corrected Transport Algorithms
for Fluids, Journal of Computational Physics, 31, pp.335‑362, 1979.
R. Leveque: High‑resolution conservative algorithms for advection in incom‑pressible flow, SIAM Journal on Numerical Analysis, 33, pp.627‑65, 1996.
[90] D. Enright, R. Fedkiw, J. Ferziger and I. Mitchell: A hybrid particle levelset method for improved interface capturing, Journal of Computational Physics,
183, pp.83‑116, 2002.
[91] http://graphics.stan ford.edu/data/3Dscanrep/

囲http://vcad‑hpsv.riken.jp/[93] H. Zhao, J. B. Freund and R. D. Moser: A fixed‑mesh method for incompressible
flow‑structure systems with丘nite solid deformations, Journal of Computational
Payszcs, 227, pp.3114‑3140, 2008.

133

囲国立天文台編:理科年表, 2008.
[95] J. C. Simo and T. J. R. Hughes: Computational Inelasticity, Springer, 1998.

[96] C. W. Macosko: Rheology principles, measurement and applications, Wiley‑
VCH Publishers, 1994.

囲R. W. Ogden: Non‑linear elastic de丘xrmations, Dover publications, 1984.

134