修士論文

シリカ充填ゴムの力学特性の有限要素均質化法による評価

指導教員：田中克志

本馬 剛徳

2015 年 2 月
神戸大学大学院 工学研究科 博士課程前期課程 機械工学専攻

Master Thesis
Characterization of Silica-Filled Rubber in Finite
Element Homogenization Method

February 2015
Division of Mechanical Engineering,
Graduate School of Engineering,
Kobe University, Kobe, Japan
TAKENORI HOMMA

要約
シリカ充填ゴムは，粒子界面の結合状態を調整する結合剤により，CB 充填ゴムより広範囲に力学特性を制御できると期待されている．本研究ではこのシリカ充填ゴムの微視的な数珠繋ぎ構造をモデル化し，ゴム相及びシリカ粒子表面のゲル相の新たな構成式の構築，均質化法に基づく有限要素シミュレーション (FEHM) によるシリカ粒子
分散構造の力学応答評価，ならびに，少数の数珠繋ぎ構造での力学特性発現メカニズ
ムの解明等を行った．
まず物理架橋点の解消を許容する非アフィン分子鎖網目理論により一般化した粘弾
性 8 鎖モデルを，化学架橋点の解消も考慮できる複合型非アフィン分子鎖網目モデルへと拡張し，さらに負荷反転時に観測される大きな応力低下を表現可能な構成式の提案も言及した．次に有限要素均質化法によりゴム相，ゲル相の絡み点数の不均一性の効果を明らかにした上で，ユニットセル中に多数のシリカ粒子をランダムに配置したモデルの等方性を明らかにした．その後，シリカ充填ゴムの巨視的応答と内部の変形に関して検討を行い，粒子に囲まれた狭い領域や引張り方向に並んだ粒子間のゴム相に大きな変形が集中している部分，粒子が凝集しているところでは粒子のクラスター
が 1 つの介在物のようにふるまうこと等を明らかにした．次にユニットセル中に 2 粒子が引張り方向に直鎖状に繋がった 0◦ ，ならびに 45◦ に傾いて繋がった 45◦ の 2 種類の代表的なモデルにより，粒子間で生じる変形集中及び粒子周りの変形メカニズムについて検討した．0◦ モデルから，粒子連結部のゲル相に非常に大きな変形が集中し，負荷変形後期における配向硬化や顕著な非アフィン変形を生じてヒステリシスロスを増大させることがわかった．一方 45◦ の場合は引張りに伴い粒子位置が中心軸方向に相対的にずれるためゲル相への変形集中は緩和され，0◦ の場合に比べユニットセル全体で非アフィン変形が発生しヒステリシスロスの増大に寄与していることがわかった．

Summary
Toward development of new FEM models that can represent the macroscopic andmicroscopic behaviors of silica-ﬁlled rubber, we proposed constitutive models for rubber and gel phase based on the non-aﬃne molecular chain network theory, and thenperformed simulations with Finite Element Homogenization Method(FEHM) with themicroscopic unit cells containing many silica particles and two-silica array.
First, we extended the previous non-aﬃne eight-chain model that accounts the eﬀectof the changes of physical/chemical entanglements during the loading and unloading
process. Additionally to reproduce the experimentally observed large stress drop at
the loading-unloading switch, a new constitutive equation accounting for the strain
rate history dependence has been proposed. Then we employ FEHM to investigate
the eﬀect of heterogeneous distribution of the entanglement in the rubber matrix andgel phase, and reestablished a simple scale-upped constitutive equation. With thus
obtained constitutive equations for the rubber matrix and gel phase, we implemented
FEHM on the microscopic unit cells with many silica particles. After clarifying theisotropy of the randomly arrayed silica-ﬁlled structure, we discussed the macroscopic
stress-strain response and internal microscopic deformation under 1 cycle loading up
to maximum stretch of λ2 =1.5. The result showed that deformation concentration
occurs at the rubber matrix in the narrow channels aggregated silica particles and also
in the wide channel between silica particles arrayed in the loading direction. We also
performed the FEHM simulation on two-silica array to discuss the fundamental aspectof deformation concentration. In the result of silica array oriented to loading axis,
extremely large stain occurs on the gel phase that connects the silica particles, andtheir non-aﬃne irreversible eﬀect contributes the large hysteresis. On the other hand,the deformation concentration was relaxed in the 45◦ inclined zig-zag silica chains.

目次
第 1 章 緒論

1

第 2 章 構成式

52.1

分子鎖網目理論 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

52.2

ゴム粘弾性体の構成式 . . . . . . . . . . . . . . . . . . . . . . . . . . .

8

2.3

非アフィン分子鎖網目モデル. . . . . . . . . . . . . . . . . . . . . . .

15

2.4

変形速度履歴依存性. . . . . . . . . . . . . . . . . . . . . . . . . . . .

18

2.4.1

変形速度履歴依存性の表現 . . . . . . . . . . . . . . . . . . . . .

18

2.4.2

未充填ゴムの変形速度履歴依存性の検討 . . . . . . . . . . . . .

20

第 3 章 粒子充填ゴムの変形応答のシミュレーション手法

23

3.1

漸近展開理論に基づく均質化手法 . . . . . . . . . . . . . . . . . . . . .

23

3.2

有限要素均質化方程式 . . . . . . . . . . . . . . . . . . . . . . . . . . .

27

第 4 章 多粒子モデルによるシリカ充填ゴムの力学特性評価

4.1

4.2

4.3

解析モデル

31

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

31

4.1.1

ゴム相及びゲル相単相の絡み点不均一性に関する均質化法解析 .

32

4.1.2

シリカ粒子分散構造の均質化法モデル. . . . . . . . . . . . . .

38

解析モデルの検証 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

41

4.2.1

モデルの等方性 . . . . . . . . . . . . . . . . . . . . . . . . . . .

41

4.2.2

解像度の検証 . . . . . . . . . . . . . . . . . . . . . . . . . . . .

42

シリカ充填ゴムの力学特性評価 . . . . . . . . . . . . . . . . . . . . . .

47

4.3.1

最大引張り時における応答 . . . . . . . . . . . . . . . . . . . . .

47

4.3.2

1 サイクルの変形過程における応答 . . . . . . . . . . . . . . . .

48

第 5 章 シリカ充填ゴムにおける変形集中機構の検討

i

51

目次

ii

5.1

解析モデル

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

51

5.2

解析結果 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

52

第 6 章 結論

60

参考文献

62

第 A 章 非圧縮性ゴム粘弾性体の構成式の速度形式表示

68

第 B 章 [ϕ],[B],[E],{ψ} の具体形

70

第 C 章 関連発表論文・講演論文

72

謝辞

77

第1章
緒論
ゴム材料は，タイヤ，衝撃吸収材，防振ゴム等の工業製品や，ボールやシューズに代表されるスポーツ用品，また医療機器等，我々の身の周りにおいて非常に幅広く用いられている．このゴム材料は，原料ゴムに数 10〜100[nm] 程度のフィラー（充填材）
を高充填することで，弾性率，引張強度，引裂き強度，破断エネルギー (1) 等の力学的特性を用途に応じて多様に変化させることが出来る．なかでも，シリカ充填ゴム（図1.1(b)）やカーボンブラック (CB) 充填ゴム（図 1.1(c)）は非常に補強効果が高いことが知られており，タイヤの材料として広範に利用されている．
ゴムのような高分子材料の内部では，モノマー（単量体）が多数連なったひも状の高分子が複雑に絡み合い三次元的な網目構造（図 1.1(d)）を形成しており，負荷時の100nm(b)Silica filled rubber.

(a)Tire.

(d)Molecular chain.

100nm(c)CB filled rubber.

Fig.1.1 Inner structure of rubber(4) .

1

応力と除荷時の応力経路が異なるヒステリシス (履歴現象)(2) や，クリープ現象，応力緩和現象，ひずみ速度依存性といった粘弾性応答 (3) 等の特異な力学特性を示す．それらの力学特性の発現メカニズムの解明を目的として，高分子材料の微視的構造の観察が進められており，柴山らは中性子小角散乱（SANS）装置を用いたゴム内部の硫黄架橋の不均一性の観察 (6) を行っている．中島らは原子間力顕微鏡（AFM）を用いた不均一構造の可視化に関する研究 (7)(8) を推進している．また池田らは変形中のゴム内部の構造変化を観察することに成功しており，分子鎖の伸長に伴う結晶化について言及している (9) ．更に分子動力学法（MD）を用いた高分子材の分子レベルの変形挙動解
析も盛んに行われており (10)(11)(12) ，屋代らは分子鎖の配向による応力上昇ならびにヒステリシス発現のメカニズムについて新たな知見を報告している (13)(14) ．
フィラー充填ゴムのヒステリシスや粘弾性応答は，前述の補強効果により未充填ゴ
ムに比べて顕著になることが確認されており (2)(15)(16) ，フィラー充填ゴムに引張変形を加えた時に生じる応力軟化の要因として，フィラー粒子の凝集体構造の変化 (17)(18)(19)や，ゴムの延伸に伴う高分子鎖の滑り (20)(21)(22) ，架橋構造の破壊 (23)(24) が古くから考えられてきた．最近では，内藤らがゴム分子鎖のエネルギー変化を詳細に検討することにより，変形中の化学架橋構造の破壊を裏付ける実験結果を報告している (25) ．
このように，フィラー充填ゴムの変形はその内部の非常に複雑なメカニズムに因るため，フィラー充填ゴムの力学特性を最大限に生かした製品設計や製造のためには，材料内部の微視的変形挙動の解明と，フィラー/ポリマー間の相互作用の適切な評価，そしてそれらを反映した巨視的応答を精密に再現し得るシミュレーションモデルの構築が必要である．我々の研究グループではこれまで，CB 充填ゴムの高機能性発現の詳細なメカニズムについて，特に CB 充填に伴うゴム部の微視的変形挙動と CB 充填ゴムの巨視的応答の関係に焦点を当てて，(1) ゴム部の適切な構成式の定式化，(2)CB 粒子が材料の機械的特性に及ぼす影響，(3) 微視領域における変形挙動を検討し得るシミュレー
ションモデルの構築及び評価など，広範多岐に及ぶ研究を推進してきた (26)(27)(28)(29) ．
そこで得られた成果はタイヤの実際の設計にも用いられている．
これまでは，工業的に利用されているゴム材料の多くは多量生産が容易な CB がフィラーとして用いられてきたが，現在は石油を原材料とせず脱石油に貢献するシリカが新たなフィラーとして期待されている．シリカ充填ゴムは CB 充填ゴムよりも耐油性・耐酸性・耐熱性に優れ，老化による架橋構造の破壊が小さいこと (30) ，またタイヤに用2

Rubber matrix.

Crosslinking Agent

Silica particle.
Gel phase.

V

V
Rubber Matrix

Rubber Matrix

50nm(a)Bunching structure.

(b)Role of crosslinking agent.

Fig.1.2 Microscopic structure of silica-ﬁlled rubber.
いる場合では転がり抵抗が小さくなるため燃費の向上をもたらし高性能なタイヤ材料となり得ることが報告されている (31) ．最近の実験では，シリカ充填ゴムの製造工程においてシリカ粒子とゴム材の界面状態を制御するために添加するシランカップリング剤（模式図 1.2(b)）の影響 (32) により，界面にゴムとは物性の異なる「ゲル相」が生成
し，図 1.2(a) のようにシリカ粒子を数珠繋ぎ状に繋ぎ止めていることが観察されている (4) ．そこで，我々の研究グループでは CB 充填ゴムを対象として開発してきたこれまでのフィラー充填ゴムの FEM モデルを基に，シリカ粒子のネットワーク構造や数珠繋ぎ構造，ゲル相の厚さなどを様々に変化させたシミュレーション (33)(34) ，また前述のゴム内部の硫黄架橋の不均一性 (6) を考慮したシミュレーション等を行ってきた (35) ．
しかしながら，実験結果に見られる引張後期の配向硬化がもたらす応力上昇，ならびに変形方向急変時における大きな応力の低下と，それに伴うヒステリシスロスが十分に表現出来ていない．更にこれらのシミュレーションで使用したユニットセルモデル
は少数の粒子を含む小規模なものであったため，実際のシリカ充填ゴム内部の構造の再現性は低く，ユニットセルの粒子配置に強く依存した基礎的な応答評価に留まっている．そこで本研究では，変形中の化学架橋構造の破壊 (25) や負荷反転時の大きな応力低下現象を表現可能な構成式を提案し，シリカ充填ゴム内部の粒子配置を模擬した多粒子解析モデルを用いた引張シミュレーションにより，フィラー/ポリマーの相互作
用，ならびにそれらを反映した巨視的応答の評価を行う．更に多粒子解析モデルでは十分な応答評価が困難であるシリカ粒子周りの解像度の良い解析モデルを用い，粒子周りの微視的変形挙動を明らかにする．
第 2 章では解析手法の基礎として，ゴムの粘弾性応答を表現するための分子鎖網目

3理論に基づく粘弾性 8 鎖モデルの構成式について説明し，変形に伴う物理架橋点及び化学架橋点の解消を表現し得る非アフィン分子鎖網目モデル，負荷から除荷へ反転する際に観察される応力低下現象を表現する構成式を提案する．第 3 章では 2 章で示したゴム材の構成式を更新ラグランジュ法に基づく有限要素均質化法として定式化する．
第 4 章では，数珠繋ぎ構造を呈するシリカ充填ゴム内部のランダムな粒子分散形態を模擬したユニットセルモデルを構築し，ユニットセルの等方性や解像度の検証を行い，シリカ充填ゴムの巨視的な変形応答の評価を行う．第 5 章では 2 粒子が直鎖状に繋がった極端な境界条件の解析モデルによる引張りシミュレーションを行い，シリカ粒子周りの局所的な変形応答について詳細に検討した．最後に，第 6 章で本研究で得られた結果の総括を述べる．

4

第2章
構成式
本章では，まずゴム弾性応答を記述するために提案された分子鎖網目理論 (36)∼(45) について説明する．次に，分子鎖に管模型 (53) を用い，非ガウス鎖理論に基づく非圧縮性を考慮したゴム粘弾性体の構成式について説明する．さらに分子鎖の絡み点数の変化を許容する非アフィンモデルへの一般化，変形方向急変時における大きな応力低下現象 (変形速度履歴依存性) を表現し得る構成式の一般化について説明する．

2.1

分子鎖網目理論

高分子とは，非常にたくさんの原子 (多くの場合は炭素原子) が共有結合によって連結したもので，図 2.1(a) に示すような長い鎖にたとえることができる．この繰り返しの構成単位をモノマーという．そして，個々の (炭素) 原子は，原子同士の結合を軸としてその周りで互いにほぼ自由に回転することができるため，全体として曲がりくねった，様々な形態をとることができる．例えば，図 2.1(a) に示す分子鎖の連続する三つの炭素原子に注目すると，共有結合による連鎖であるから，結合長さ l = 1.54Å，結合角 θ = 70.53◦ と確定している (36) ．これに対して，第 4 番目の炭素原子の結合は，l と

θ を一定に保ちながら，第 2 結合を軸に回転可能となり，その位置は，回転角の関数として表されるポテンシャルエネルギによって決まる．このような考え方で高分子材料の微視的構造を忠実に考慮したモデルを構築し，高分子材料の挙動を表現することが原理的には可能である．しかしながら，実際の適用に当たっては多くの時空的な制約が加わるため，モノマーを直接扱わず，セグメントという最小構成単位で粗視化した分子鎖網目モデル (37) が用いられる．

5

Fig.2.1 Concept of hierarchical structure of polymeric material(37) ,(a)molecular chain, (b)segment, (c)chain network structure,(d)macroscopic continuum.

分子鎖網目理論では，高分子材料は図 2.1(c) に模式的に示すように，分子間の化学的結合あるいは物理的結合により接合点において連結された鎖が，ランダムに配向した網目構造を有していると仮定している．さらに，(i) 接合点は原子の揺らぎ周期に対して長時間的には平均位置が変化せず，接合点周りの摂動は無視できる，(ii) 二つの接合点を両端に持つ分子鎖の端−端ベクトル (end-to-end vector) は，それが埋め込まれている材料の連続体と共変形をするとの仮定を置く．このようなモデルをアフィンモ
デルという．図 2.1(b) に示すように，二つの接合点間の分子鎖は「1 本の分子鎖」と定義され，それは導入した最小構成単位である「セグメント」から成る．モノマーの数が十分多ければスケーリング則によって鎖の巨視的な性質は変わらない (39) ．ここで導入した「セグメント」が「モノマー」とどう対応するかは議論の対象とせず，ここでは現象論的に最小単位のセグメントを導入している．
「1 本の分子鎖」の形態が非ガ
ウス統計分布 (40) に従うとすると，二つの接合点を結ぶ方向にストレッチ λ を加えた場合に生じる応力 σ は次式で表すことができる (38) ．
)
(
√
λ
−1
√
σ = kB T N λL
N

(2.1)

ここで，N は 1 分子鎖あたりのセグメント数，kB は Boltzmann 定数，T は絶対温度で6

Fig.2.2 Molecular chain network model, (a) three chain model, (b) four chain
model, (c) eight chain model.
ある．また，関数 L(x) は次式で定義される Langevin 関数である．
{ ()}
sinh xd
1
L(x) =ln= coth x −
dx
x
x

(2.2)

網目の全体的な応答特性は，個々の鎖の寄与を考えることにより得ることができるが，その取り扱いは数学的に極めて困難なものとなる．そこで，網目構造の応答モデルを得るために簡便な平均化手法が提案されている．

James 及び Guth(41) は単位体積あたり n 本の鎖を含む網目は直行する 3 本の軸方向に n/3 ずつの鎖が配置されたものと相当であると仮定した，図 2.2(a) に示す 3 鎖モデルを提案した．Wang 及び Guth(42) はこの 3 鎖モデルを等 2 軸変形に適用した．同様に，Treloar(43) は，図 2.2(b) に示す 4 鎖網目モデルの概念 (44) をゴム弾性に適用することを提案したが，主ひずみ空間における対称性を表現することができないことが示されている (45) ．Arruda 及び Boyce(45) は図 2.2(c) に示す 8 鎖モデルを提案し，これらの網目モデルの中で最も広範な変形モードに適用できることを示した．

3 鎖モデルの場合，各主ストレッチ方向の変形は互いに干渉しないため主応力はそれぞれの軸方向での応力の足し合わせにより計算できるが，8 鎖モデルの場合は 3 軸方向の変形がそれぞれ干渉し合うため，1 本の分子鎖方向のストレッチ λc は，主スト√
レッチを λ1 ，λ2 ，λ3 とすると，λc = (λ21 + λ22 + λ23 ) /3 と表すことができるので，

∂λc
λi
=
∂λi
3λc
の関係が成り立つ．

7

(2.3)

一方，変形前の体積を基準にした単位体積あたりの仕事に相当するひずみエネルギー密度関数 W を用いると，1 分子鎖の応力 σc は次式のように表せる (48) ．

σ c = λc

∂W
∂λc

(2.4)

式 (2.1) のストレッチ λ を分子鎖方向のストレッチ λc に置き換え，さらに式 (2.4) より恒等的に次式が成り立つ．

√∂Wc= kB T N λc L−1
∂λc

(

λ
√c
N

)
(2.5)

式 (2.3)，式 (2.4)，式 (2.5) より，8 鎖モデルの主ストレッチ方向の応力 σi とストレッチ λi は次の関係で与えられる．

∂W∂W ∂λc= λi∂λ∂λc ∂λi
{ i
(
)} 2
√
λc
λi
1
−1
√
=
kB T N L
3λc
Nσi = λi(2.6)

本研究では，この 8 鎖モデルを基礎としてゴム粘弾性体の構成式を定式化する．

2.2

ゴム粘弾性体の構成式

図 2.2(c) で示される８鎖モデルはゴム超弾性体の変形応答を記述するのに用いられている．超弾性体とは負荷を受け大きく変形した後，完全に除荷すると元の状態に戻る弾性体である (46) ．しかしながら，実際の分子鎖は周囲の分子鎖からの摩擦に起因する粘性も持ち合わせている．そこで，ゴムの粘弾性挙動を記述するために，図 2.3 に示すような粘弾性８鎖モデルとダンパーで構成されるモデルを構築する．
周囲の分子鎖との相互作用を表現するために，図 2.2(c) に示す 8 鎖モデルの各単鎖に，図 2.3 挿入図に示す粘性抵抗をもつバネ・ダンパーの標準モデルを導入した新たな 8 鎖モデル A が提案されている (47) ．ここで用いた粘性抵抗は，後に説明する管模型 (53) によって表される．本モデルにおける応力の評価はバネ部分のみで行い，バネとダンパーを直列関係にすることで粘性効果を応力の評価に反映させている．ここで，初めに８鎖モデル A の構成式を記述する．8 鎖モデルの要素 α，β ，γ のストレッチをそれぞれ λα ，λβ ，λγ とすると，1 本の分子鎖の変形量と要素 α の変形量は等しく，ま

8

B
A

β
α
γ

D

Fig.2.3 New visco-elastic model consists of revised 8 chain model A, conventional eight chain model B and dashpod D.
た要素 β ，γ は要素 α と並列の関係にあることから，以下の関係が成り立つ．

λc = λα = λβ λγ(2.7)

式 (2.6) を 8 鎖モデル A に適用すると，8 鎖モデル A の主ストレッチ方向の応力 σiA とストレッチ λi は次の関係で与えられる．

1
σiA =
3

(
)}
{ √
(
)
√
λ
1
λ
λ2i
c
β
CαR Nα L−1 √+CβR Nβ L−1 √
λγ
λc
Nα
Nβ

ここで CαR = nα kB T ，CβR = nβ kB T ，Nα ，Nβ は分子鎖のセグメント数で

(2.8)

√
√
Nα ， Nβは分子鎖の限界伸びを表す．nα ，nβ は構成要素のバネ α，β に含まれる分子鎖の数を表す．添え字 α，β ，γ については，ストレッチと同様図 2.3 に示す要素 α，β ，γ と対応している．
ここで弾性項のみである図 2.3 の 8 鎖モデル B の主ストレッチ方向の応力 σiB と弾性ストレッチ λ′i は次のように表せる (45) .
{()} 2 ′
√
1
λcB
λi
B
−1
R
√
σi =
CαB NαB L
3
λcB
NαB

(2.9)

R
CαB
= nB kB T ，n = nα + nβ + nB ，n は単位体積中に含まれる鎖の数を表す．一般に，

ゴム粘弾性体の変形は体積変化が小さいとしてそれを無視する場合が多い．そこで本
研究では，非圧縮性ゴム粘弾性体を取り扱うものとし，非圧縮性を満たすために静水

9圧 p を用いる．この時，式 (2.8) , (2.9) を用いると，非圧縮性ゴム粘弾性体の構成式は次式のように表せる．

σi = σiA + σiB − p(2.10)∇

また，構成式 (2.10) の速度形式は，Kirchhoﬀ 応力の Jaumann 速度 S ij とひずみ速度テ
ンソル ε̇kl ，粘性ひずみ速度テンソル ε̇pkl を用いて，次のように表すことができる．

[{
(
)}
√
)
R
′
′
√ ( ζ
C
N
1
L
L
ζ
β
β
√
CαR Nα √−+−Aij Akl /Amm
S ij =3λγNα λc
λγ Nβ λc]{√ }√LCαR Nα L′ CβR Nβ++{δik Ajl + Aik δjl } ε̇kl
λc
λc
(
)
√
[{
)}
( ′′
′
√
CβR Nβ λ̇γ
λ
ζ
ζ
L′′
1
β
′
R
L +√
CαB NαB √−A′ij A′kl /A′mm
− 2√
Aij +
3
λγ 3AmmNαB λcBNβ√]R}
L′′ CαBNαB {′′+δik Ajl + Aik δjl (ε̇kl − ε̇pkl ) − ṗδij(2.11)λcB∇

式 (2.11) の具体的な導出方法については [付録Ａ] を参照されたい．ここで，Aij は√√
左 Cauchy-Green 変形テンソル，L = L−1 (λc / Nα )，L′ = L−1 (λβ / Nβ )，L′′ =√L−1 (λcB / NαB )，ζ = L2 /(1 − L2 csch2 L)，ζ ′ = L′ 2 /(1 − L′ 2 csch2 L′ )，ζ ′′ = L′′ 2 /(1 −L′′ 2 csch2 L′′ ) である．また，添え字に図 2.3 と対応したものを付した．本研究ではペナ
ルティ法を用いることにより，非圧縮性を近似的に満足させる．
次に粘性抵抗を表現する，要素 β ，γ の扱いについて説明する．ここでは，要素 B，

D も要素 β ，γ と同様の動きをすると仮定して取り扱う．ある時刻に変形勾配が F となるような負荷あるいは変形を受けているゴムの変形を考える．その時の変形勾配 Fを次式で定義する (49) ．

F =

∂x
∂X

(2.12)

ここで，X は物体点の基準配置，x は現在の配置を表す．ゴムにおける基準配置は分子鎖がランダムに配向した等方性状態である．図 2.4 に示すように，変形勾配 F は弾性部分 F β と粘性部分 F γ に次式のように分解できる．

F = F βF γ

(2.13)F γ は完全な除荷状態で応力解放配置を表す．また，変形勾配 F は弾性ストレッチ V β ，10FCurrent
Configuration
Reference
Configuration

F

Fβ

β-1

RelaxedConfiguration
γ

F

Fig.2.4 Concept of viscoelasitc decomposition of deformation gradient.

回転 R，粘性ストレッチ U γ を用いて次の形で表現される．

F = V β RU γ

(2.14)

回転 R を弾性部分と粘性部分に分け，

R = Rβ Rγ ,(2.15)

極分解定理 (50, 51) に従うと，次式の関係を得る．

F β = V β Rβ = Rβ U β(2.16)

F γ = Rγ U γ = V γ Rγ

(2.17)

実際，回転は弾性か粘性かは特定することはできない．しかしながら，ここではRβ = I,R = Rγ(2.18)

F γ = V γ R = RU γ

(2.19)とすることによって，次式の関係を得る．

11

つぎに，速度勾配 L を考える．

L=∂v= d + w = Ḟ F −1 = Ḟ β F β−1 + F β Ḟ γ F γ−1 F β−1∂x(2.20)

ここで，v は変位速度，d は変形速度テンソルで L の対称部分，w はスピンテンソル
で L の反対称部分である．また，d と w をそれぞれ弾性成分と粘性成分の和であるとすると，次の表現が得られる．

d = dβ + dγ ,dβ + wβ = Ḟ β F β−1 ,w = wβ + wγ
dγ + wγ = F β Ḟ γ F γ−1 F β−1(2.21)(2.22)

応力解放配置の速度勾配 Lp は，次式で与えられる．
γ

Lγ = Ḟ γ F γ−1 = d̃ + w̃γ(2.23)式 (2.1) より，要素 β の応力 σβ とストレッチ λβ のは，次の関係で与えられる．
(
)
√
λ
β
(2.24)σβ = CβR Nβ λβ L−1 √Nβ
ただし，σβ = σγ(2.25)σ γ は要素 γ にかかる応力である．式 (2.7) より，弾性ストレッチ λβ は分子鎖一本のス
トレッチ λc を用いて，次のように表される．

λβ =

λc
λγ

(2.26)

粘性ストレッチ λγ は次のように表現される．

∫ t
λ̇γ dtλγ =(2.27)0˜ γ|
λ̇γ |t=t+∆t = λγ |t=t ·|d(2.28)

粘性変形速度 d̃γ は負荷あるいは除荷配置のどちらの場合においても一般的に次のように表現されるべきであると考える (49, 52)d̃γ = γ̇ γ N
12(2.29)

ここで，γ̇ γ は粘性せん断ひずみ速度，N は方向を示すテンソルである．粘性流れの駆
動応力 (driving stress) σ ∗γ は Cauchy の応力 σ γ を用いて次のように表される．

σ ∗γ = σ γ(2.30)

連合流れ則によって粘性変形速度が偏差駆動応力方向に発生すると仮定すると N は次式のようになる．

σ ∗′
γ
N=√
2 τ∗

(2.31)

ここで，( )′ は偏差成分を表し，τ ∗ は N を単位の値として定義するために導入した量で，次のように表すことができる．

[
∗

τ =

1 ∗′ ∗′
σ · σγ2 γ] 12(2.32)

Doi·Edwards(53) は，高分子鎖で起こる特異な粘弾性現象を説明するために，de Gennes(54)によって提案された管模型 (tube model) を分子鎖に適用した reptation (爬行) 理論を示した．図 2.5 に示す管模型では，周囲の分子鎖との摩擦を，分子鎖の主鎖と直交方向の運動の制限と捉え，分子鎖の主鎖方向への運動は自由であるがその垂直方向への運動は周囲の分子鎖にあまり影響を与えない程度の距離 a，長さ L の管内で拘束されていると仮定している．外力を加えると管は変形し，管の直径方向，軸方向ともにまず a の距離内にあるセグメントの配向分布の緩和が短時間のうちに起こる．一方，軸方向には reptation 運動により分子鎖は最初のゆがんだ形状の管から徐々に抜け出し，aaFig.2.5 Tube model.

13Undeformed Network

Deformed Network

Deformed and Relaxed Network

Fig.2.6 Relaxation behavior of polymer chain.

完全に抜け出したとき，すなわち分子鎖が管に沿って長さ L だけ移動したとき応力は完全に緩和する．このような周囲の分子鎖との摩擦によって，実際の分子鎖網目構造
では図 2.6 に示すような緩和現象が生じていると考えられる．

Bergström·Boyce は reptation 理論を基に，与えられた有効せん断応力 τ ∗ に対する√
粘性せん断ひずみ速度 γ̇ γ を次のように導出した (3) ．分子鎖の変位を û = a1 ϕ(t) とする時，緩和時間 t と分子鎖長 l(t) の関係は

l(t) = l0 + a1√
ϕ(t)(2.33)と表される．ここで l0 は初期分子鎖長，ϕ(t) は reptation 理論による緩和時間である．
粘性ストレッチ λγ は次式のように表される．

λγ (t) =l(t)= 1 + a2 ta3
l0(2.34)

ここで a2 > 0，0.5 < a3 < 1.0 である．時間微分を取ると，次のようになる．

λ̇γ = a2 a3 ta3 −1(2.35)式 (2.34)，(2.35) よりクリープ速度は次式のように表される．

λ̇γ = a4 (λγ − 1)a5

(2.36)

ここで，a4 > 0，a5 ∼= −1 である．しかし，クリープ速度は駆動応力に依存するとされるので，粘性せん断ひずみ速度を以下のように表した．

γ̇ γ = Ĉ1 [λγ − 1]C2 τ ∗ m
14(2.37)

ここで，Ĉ1 ，C2 ，m は材料定数であり，一般にひずみ速度に依存する．

2.3

非アフィン分子鎖網目モデル

ゴムは高分子鎖がランダムに結合した網目構造を有し，網目の接合点として振る舞
う絡み点は分子間の共有結合による化学架橋点と，それに比べ結合力の弱い分子間力によって結合している物理架橋点に分類できる．系の中で絡み点数が多いということ
は，絡み点間の分子鎖長さが短い，すなわち 1 分子鎖当たりの平均セグメント数 N が小さいことに対応する．変形過程において分子鎖が滑り出すと架橋点が消滅し，絡み点数が変化することが実験的に示唆されており (55),(56) ，対応した非アフィン分子鎖網目モデルが提案されている (57),(58) ．網目構造の変形で絡み点数が減少することによって，1 分子鎖あたりの平均セグメント数 N は増加する．絡み点間の分子鎖を「1 本」とカウントしているため，図 2.7 に模式的に示すように絡み点が解消すると単位体積中の分子鎖数 n(= nα + nβ + nB ) は減少し，伸長可能性の向上と剛性の低下をもたらす．
本研究では，セグメント数 N が分子鎖ストレッチ λc に依存すると仮定して，Junction point

1

Junction point21(a)Before(b)AfterFig.2.7 Schematic view of extinction of junction point by decrease in physicallinkage.

N (λc ) = N0 + f (λc )(2.38)

Nｎ = nα Nα + nβ Nβ + nB NB = constant(2.39)

15

CαR
nα =
,
kB T

CβR
nβ =
,
kB T

R
CαB
nB =
kB T

(2.40)とする．ここで，Nn は系全体の総セグメント数，N0 は初期セグメント数であり，総セグメント数 Nn は材料に固有の数値である．
次に関数 f (λc ) について説明する．一般に，高分子材料の変形過程において消滅す
るのは結合力の弱い物理架橋点のみであると考えられており，現在提案されている非アフィン分子鎖網目モデルは (57),(58) ，この物理架橋点のみの消滅に対応している．しかしながら，最近の研究報告により高分子材料の変形過程の初期段階において，化学
架橋点の消滅が顕著になる可能性が示唆されている (25) ．そこで本研究では物理架橋
点，化学架橋点それぞれの消滅に対応した関数 f (λc ) の構築を行い，それらを混合して表現出来るように拡張した複合型非アフィン分子鎖網目モデルを提案する．
まず，物理架橋点の消滅に伴う 1 分子鎖当たりのセグメント数 N の変化を表す関数

fph (λc ) を下記のように，
fph (λc ) = a0 + a1 λc + a2 λ2c(2.41)とする．a0 ，a1 ，a2 は定数である．次に化学架橋点の消滅は変形初期に多く発生することが示唆されていることから，関数 fch (λc ) を下記のように，

fch (λc ) = a0 + a1 log(λc − 1)

(2.42)と対数を用いて表現することとした．更に物理架橋点と化学架橋点の消滅を混合して表すために，それぞれの影響を重み関数 w で制御する次の複合型非アフィンモデルを構築した．

fmix (λc ) = wfph (λc ) + (1 − w)fch (λc )(2.43)

w の値が 0 に近づくほど化学架橋点消滅の影響が大きくなり，反対に 1 に近づくほど物理架橋点消滅の影響が大きくなる．図 2.8(a) に重み関数の値を様々に変化させた
時の平均セグメント数 N とストレッチ λ の関係を示す．こちらに示すのは最大引張り時のストレッチを 4 とし，重み関数 w の値に因らず初期セグメント数 N 及びストレッチが 3 の時の平均セグメント数 N を統一した時の 1 サイクル分の N の変化である．負16

True Stress σt22[MPa]Number of segment N

25

20

15

w=0.0
w=0.25w=0.5
w=0.75
w=1.0

10
1

2

3

6

4

2

0

4

w=0.0
w=0.25w=0.5
w=0.75
w=1.0

1

2

3

4

Stretch λ(b)Stretch λ(a)Fig.2.8 Comparison of (a)number of segment - stretch relations and (b)true
stress - stretch relations by weighting function.

荷から除荷に移行してからの絡み点数変化は不可逆なものとし，本シミュレーション
では 2 サイクル目以降の再負荷時において，前のサイクルで受けた最大ストレッチに達しない変形領域では平均セグメント数 N は変化しないものとした．更に変形が進み，前のサイクルで受けた最大ストレッチを超えると，再び式 (2.38) により N が変化する
ものとする．また図 2.8(b) には前節で示した 8 鎖モデル A の計算式 (2.8) の Nα にセグ
メント数変化式 (2.38) を導入して計算した結果の真応力-ストレッチ関係を，それぞれ
図 2.8(a) のセグメント数変化に対応させて示す．構成式による応力計算では応力値はセグメント数 N とストレッチ λ にのみ依存するため，セグメント数 N の一致するス
トレッチが 3 の時点での応力値は一致しており，この点を境に応力値の大小が入れ替わっている．重み関数 w の値による特徴として，w が 0 に近いモデルでは負荷変形後期においても応力は軟化せず，1 サイクルの変形におけるヒステリシスロスは小さい．
反対に w が 1 に近いモデルでは負荷変形後期において応力軟化が見られ，ヒステリシ
スロスが大きくなる傾向にある．本研究ではこれらの特徴を踏まえた上で，未充填ゴ
ムの実験結果の負荷時における負荷変形中の応力軟化や負荷変形後期の配向硬化，ヒステリシスロスを再現出来るパラメータを決定する．

17

2.4

変形速度履歴依存性

2.4.1

変形速度履歴依存性の表現

ゴム等の高分子材料の変形過程におけるヒステリシス発現のメカニズムとしては，

2.2 節に示した粘弾性現象や 2.3 節に示した非アフィン変形が大きな要因として考えられている (20)(23) ．しかしながら粘弾性現象が顕在化しないほどの低ひずみ速度や非ア
フィン変形による絡み点解消の影響が小さい低ひずみ領域においてもヒステリシスの発現が見られることから，ヒステリシスの発現には上記に加えて別の要因が存在する
ことが示唆される．そこで本研究では，ゴムの繰り返し変形過程において，負荷方向が急変する場合大きな応力変化が発生し，それが変形の進行とともに緩和されていく
実験事実に注目した．この現象を，変形方向が急変するとゴム内部の 1 本の分子鎖に対する周囲の分子鎖からの拘束が強くなり，変形方向急変直後に少ない変形量で大きな応力の変化が生じる，すなわち分子鎖の剛性が一時的に上昇したとみなして表現する．
これを「変形速度履歴依存性」として 2.2 節のゴム粘弾性体の構成式モデル (図 2.3) の弾性変形を担う α の部分，すなわち式 (2.8) の第 1 項を変形速度履歴依存性が表現出来るように一般化した．
本研究で用いる構成式 (2.8) 第 1 項においてゴムの剛性を支配するパラメータは CαR ，√
Nα の 2 つであり，どちらも 2.3 節で説明した非アフィン変形により変化するパラメー
タである．そこで本研究ではこれらのパラメータに加え，新たなゴムの剛性パラメー
タとして変形速度履歴パラメータ ξ を導入し，式 (2.44) に示すように構成式を変形速度履歴依存性が表現出来る形式に一般化した．

1
σiA =
3

{

)
(
√
√
1
λc
R
−1
√
Cα Nα (1 + ξ)L+CβR Nβ L−1
λγ
Nα

(

λ
√β
Nβ

)}

λ2i
(2.44)λc

ξ は負荷時には 0 とし，また変形速度履歴依存性は低ひずみ速度条件下，つまりひずみ速度には無依存であるため分子鎖ストレッチ λc にのみ依存するものとして次式で表す．

ξ = a0 + a1 λc + a2 λ2c
18(2.45)

1

(3)
0.5

ξ

(2)

(1)
0

1

2

3

4

Stretch λFig.2.9 ξ - stretch relations.

1 サイクルの変形過程における ξ の値の変化とストレッチ λ の関係の一例を図 2.9 に示す．ここで，図の横軸は λc ではなく λ であることに注意されたい．図中の矢印で示すように変形過程において (1) 負荷時には ξ の値は 0 のまま推移し，(2) 負荷反転時に大きな ξ の値をとることでゴムの剛性の上昇を現象論的に表現する．(3) そしてその後除荷が進むにつれ剛性が減少し，ξ の値が一時的に負になるようにすることで実験結果を表現することを試みた．
ただし，変形速度履歴依存性がゴムの変形応答に与える影響の大きさについては実験的に明らかにされていないため，変形方向急変時の ξ のピーク値は 1，あるいは 2(剛性はもとの 2 倍，あるいは 3 倍) と仮定する．本節ではこの ξ を導入した影響が低ひず
み速度，低ひずみ領域においても発現することを確認するために，ひずみ速度による比較，最大引張り時のストレッチによる比較，ξ の値による比較をそれぞれ式 (2.44) を用いて行う．

2.4.2

未充填ゴムの変形速度履歴依存性の検討

図 2.10 にひずみ速度 ϵ̇ = 10−5 ，ϵ̇ = 10−1 の時の 1 サイクルの変形における (a) 真応
力-ストレッチ関係，(b) ヒステリシスロスを示す．なお，ひずみ速度以外の計算条件は以下全ての計算において表 2.1 に示すものとし，重み関数 w は 0.5 とする．ξ の値はひずみ速度には無依存であるため，ひずみ速度の大小に関わらず図 2.9 に示すように推移

19

Table 2.1 Parameter of revised new visco-elastic model.

CαR (= nα kB T )0.275

CβR = (nβ kB T )0.40

CBR (= nB kB T )
0.50

Nα
12.0

Nβ
12.0

NB
12.0

C1γ
3.0 × 102

C2γ
−0.50

mγ
6.4

C1D
3.0 × 103

C2D
−0.50

mD
4.8

していくものとしている．ここでまず ξ の値は除荷時にしか変化しないため負荷時の応力曲線は変形速度履歴無依存の場合と同じ経路を辿っている．次に負荷反転後，除荷時の応力曲線については ξ を導入したモデルの方が大きく応力が低下しており，それに伴いヒステリシスロスが増大していることがわかる．またヒステリシスロスの増加率を比較すると，ひずみ速度 ϵ̇ = 10−5 の場合は約 1.88 倍，ϵ̇ = 10−1 の場合は約 1.44倍となっており，ひずみ速度が小さい時の方が ξ の影響が大きくなっている．これは式 (2.44) の 2 つの項のうち ξ を導入している第 1 項，つまり応力全体に対する弾性変形部分の担う応力が大きいためである．以上のことから ξ を導入することで低ひずみ速度の条件下であっても負荷反転時の大きな応力低下やヒステリシスロスの増大を表現出来ることが確認出来た．しかしながら，除荷終了時に着目すると非常に大きな残留ひずみが発生していることがわかる．ゴムなどの高分子材料の繰り返し変形において，ひずみ速度が大きいと粘性の影響による応答遅れで除荷終了時には残留ひずみが発生することが知られているが，ξ を導入すると低ひずみ速度の場合でも残留ひずみが非常に大きくなってしまうことがわかった．
次に図 2.11(a) に最大引張り時のストレッチを λ = 2.0, λ = 3.0, λ = 4.0 の 3 パターンに変化させた時の ξ とストレッチの関係を示す．計算条件としてひずみ速度は ϵ̇ = 10−1とし，変形速度履歴依存性を低ひずみ領域でも発現することから図 2.11(a) に示すように最大引張り時のストレッチの大きさに関わらず ξ の値のピークを 1 で統一している．
また図 2.11(b) に真応力-ストレッチ関係を示す．最大引張り時のストレッチが λ = 4.0の高ひずみ領域においては除荷時の大きな応力低下が見られるため ξ の影響が確認出
来るが，低ひずみ領域においては ξ の有無による大きな応力低下は表現出来ていない．
このことから ξ を導入した場合でも低ひずみ領域における変形速度履歴依存性，すな
わち大きな応力の低下の表現が不十分であることがわかった．

20

10Previous model(ε = 10−5)
Invert loading parameter model(ε = 10−5)
Previous model(ε = 10−1)
Invert loading parameter model(ε = 10−1)

8Previous model(έ=10-5)
Invert loading parameter model(έ=10-5)

Hysteresis loss [J/m3]

True Stress σt22[MPa]

10

6

4

2

0

1

2

3

8

Previous model(έ=10-1)
Invert loading parameter model(έ=10-1)

6

4

2

0

4

Stretch λ(b)(a)

Fig.2.10 Comparison of true stress - stretch relations by strain rate.

102

λmax=2.0λmax=3.0
λmax=4.0

True Stress σt22[MPa]Previous model
Invert loading parameter model

ξ

1

0

1

2

3

8

6

4

2

0

4

1

2

3

Stretch λStretch λ

(a)

(b)

4

Fig.2.11 Comparison of (a) parameter ξ-stretch relations and (b) true stress stretch relations by max stretch.

次に図 2.12(a) に ξ の値のピークを 0(=変形速度履歴無依存)，1，2 と変化させた時の ξ とストレッチの関係を，(b) に真応力-ストレッチ関係を示す．ξ の値を大きくすると除荷時の大きな応力の低下は非常に顕著になり，それに伴いヒステリシスロスも増

21

10

2

True Stress σt22[MPa]ξmax=0.0
ξmax=1.0
ξmax=2.0

ξ

1

0

1

2

3

8

6

4

2

0

4

Previous model
Invert loading parameter model(ξmax=1.0)
Invert loading parameter model(ξmax=2.0)

1

2

3

Stretch λStretch λ

(a)

(b)

4

Fig.2.12 Comparison of (a) parameter ξ - stretch relations and (b) true stress- stretch relations.

大するが，ひずみ速度を変化させた場合と同様，残留ひずみが大きくなりすぎることがわかった．
以上の検討から，変形速度履歴依存性の導入によるひずみ速度の大小に関わらず負荷反転時の大きな応力の低下やヒステリシスロスの表現は良好であったが，低ひずみ領域では不十分であることや除荷終了時の残留ひずみの過大評価等の問題点が浮上した．本論文の 4 章以降での有限要素均質化法によるシリカ充填ゴムの解析では，ストレッチが λ = 2.0 程度の低ひずみ領域で変形している要素が大部分を占めるため，ξ の
影響はほとんどない．したがって，変形速度履歴依存性を表現するための新たな構成
式の提案は本章で終わりとし，4 章以降の解析には変形速度履歴依存性を考慮しない
構成式を用いて進める．

22

第3章
粒子充填ゴムの変形応答のシミュ
レーション手法
本章では，2 変数漸近展開理論に基づく均質化法の基本的な考え方を述べた後，ゴムの構成式を更新ラグランジュ法に基づく均質化理論 (59, 60) に導入することにより，微視的関係式及び巨視的平衡式を導出し，その有限要素表示式を示す．

3.1

漸近展開理論に基づく均質化手法

本節では，2 変数漸近展開理論に基づく均質化法の基本的な考え方を簡単に述べた
のち，ゴム材の構成式を更新ラグランジュ法に基づく均質化法 (59, 60) に導入することにより，微視的関係式及び巨視的平衡式を導出する．
図 3.1 に示すような全体構造 X の任意点の近傍において，局所的に周期性をもつ

W
xi

P

X

unit cellStYSu

yi = xi /hyi

Fig.3.1 The relation between macroscopic continuum body and internal microscopic structure.

23

微視構造 Y が存在する材料を仮定し，構造物全体を表現する座標系 xi (i = 1, 2, 3) とyi = xi /η の関係を満足する微視構造を表現する座標系 yi の 2 変数を導入する．ここで，η は微視的周期構造内の基本単位領域のスケールを表す．現変形状態における物体の体積を Ω，表面積を S ，外部表面の一部 St 上に作用する表面力を P とし，残りの外部表面 Su に一定の変位速度を与える．このとき，更新ラグランジュ法を用いると，
仮想仕事原理式は下記のように表せる (61) ．

∫ (

∫

)
Ṡji + σmj vi,m δvi,j dV =VṖi δvi dS

(3.1)

St

ただし，Ṡij は Kirchhoﬀ の応力速度を表す．一方，式 (2.11) の構成式は下記の関係式(61)∇

Ṡij = S ij −Fijkl ε̇kl ,1
Fijkl = (σlj δki + σkj δli + σli δkj + σki δlj )
2

(3.2)

を用いることによって，次のような形に統一的に示すことができる．

Ṡij = Lijkl ε̇kl

(3.3)

なお，Lijkl は 2 章で示した速度勾配である．一方，微視領域内の任意点の変位速度v はスケールパラメータ η により，次のように漸近展開できる．

v = v (x, y) = v 0 (x, y) + ηv 1 (x, y) + η 2 v 2 (x, y) · · ·(3.4)

式 (3.4) をひずみ速度と変位速度の関係式

1ε̇ij =
2

(

∂vj
∂vi
+
∂xj ∂xi)(3.5)に代入し，次式が得られる．

][
1ε̇ij = ė0ij (v) + Ėij0 (v) ė1ij (v) + η Ėij1 (v) + ė2ij (v) · · · ,
η
(
)
(
)
k
k
k
k
∂v
∂v
1
∂v
1
∂v
j
j
i
i
ėkij (v) =++, Ėijk (v) =
2 ∂yj
∂yi
2 ∂xj
∂xi
24

(3.6)

次に，式 (3.3)，(3.4)，(3.6) を式 (3.1) に代入し，η について同じ次数の項を整理すると，以下の式が得られる．

)
∫ (
∂δvi
∂vi0 ∂δvi
0
Lijkl ėij (v)+ σmj
dV = 0∂yj∂ym ∂yj
Ω

1
η2

1
η

∫ {

[(

) ∂δv

0
(v) + ė1kl (v)
Ėkl

∂δvi
+ ė0kl (v)
∂yj
∂xj
]}
0
∂vi ∂δvi

∫ {
Ω

[(
Lijkl
[(

+σmj

) ∂δv

0
Ėkl
(v) + ė1kl (v)

(

i

]

i

Lijkl
[( 0
)
∂vi1 ∂δvi
∂vi
+σmj
+
+
∂xm ∂ym ∂yj∂ym ∂xjΩ(3.7)

dV = 0(3.8)) ∂δv ]

1Ėkl(v) + ė2kl (v)

i

+
(3.9)
∂xj
∂yj
)
( 1
)
]}
∫
∂vi
∂vi0
∂vi1 ∂δvi∂vi2 ∂δvi
+
+
+
dV =
Ṗi δvi dS∂xm ∂ym ∂xj∂xm ∂ym ∂yj
St

一方，Y-periodic 条件を満たす関数 Ψ (y) に対して，

( )
∫ ∫
x
1
Ψ (y) dY dΩ,ΨdΩ →η|Y | Ω YΩ(3.10)( )∫ ∫x1lim η ΨdS →Ψ (y) dSdΩη→0+η|Y | Ω S
S(3.11)

∫
lim

η→0+

∫

が成立する (62) ．ここで，|Y | は微視領域の体積である．式 (3.10) を用い，式 (3.7) から次式が得られる．

()]∫ {∫ [
∂
∂vk0
∂vi0
1
−
Lijkl
+ σmj
δvi dY|Y | Ω
∂yj
∂yl
∂ym
Y
)
}
∫ (
∂vk0
∂vi0
+
Lijkl
+ σmj
nj δvi dS dΩ = 0
∂yl
∂ym
S

(3.12)

δvi が任意であるため，式 (3.12) から次式が得られる．
∂
−
∂yj
(

(
)
∂vk0
∂vi0
Lijkl
+ σmj
=0
∂yl
∂ym

∂v 0
∂v 0
Lijkl k + σmj i
∂yl
∂ym
25

(3.13))
nj = 0(3.14)

Guedes ら (62) の命題 1 に基づき，式 (3.13)，(3.14) から，次式を得る．

v 0 = v 0 (x)(3.15)

式 (3.15) より，変位速度の漸近展開式 (3.4) の第一項 v 0 は巨視的な座標系 x にのみ依存することが分かる．次に，式 (3.8) に対し，式 (3.10)，(3.11) を用いると，次式を得る．

∫∂v 1 ∂δvi
dY = −(Lijkl + σlj δik ) k∂yl ∂yjY∫(Lijkl + σlj δik )Y∂vk0 ∂δvi
dY∂xl ∂yj(3.16)0

式 (3.16) は v 0 に対して線形であるので，v 1 と Ė は次式に示す関係が存在する (62) ．
0v 1 = χĖ (v)(3.17)

ただし，χ は特性変位関数と呼ばれる Y-periodic を満足する関数で，それぞれは下記の式の解である．

(
)
]
∫ [
∫
kl
∂χ
1 ∂χkl
∂χkl
∂δvi
∂δvi
p
p
m
Lijpm
+
+ σmj δpidY =(Lijkl + σlj δki )
dY
2 ∂ym
∂yp
∂ym ∂yj∂yjyY(3.18)
さらに，式 (3.10) において，可容変位速度 δv は任意に選ぶことができるので，δv =δv (x) とし，式 (3.17) を用いることにより，次式が得られる．
∫ [

]

∫
∂δvi
dΩ =
Ṗi δvi dS,
∂xj
Ω
St
)]
(
∫ [
kl
kl
∂χ
∂χ
1
1
p
q
Lijkl − Lijpq+
dY,
LH
ijkl =|Y | Y
2 ∂yq
∂yp
)
∫ (
1
∂χkl
i
H
τijkl =dYσlj δki − σmj|Y | Y
∂ym
0
0
H ∂vk
Ė
(v)
+
τ
LH
ijkl kl
ijkl∂xl(3.19)

以上から，微視構造について解くべき特性変位関数 χ は，微視構造の形態と材料定数のみに依存し，全体構造のひずみ，応力などから独立して求解されることが分かる．
一方，全体構造について解くべき巨視的平衡方程式 (3.19) は均質化された巨視的特性
量などが特性変位関数より求められるため，微視構造と独立して求解することが可能
となる．

26

有限要素均質化方程式

3.2

本節では，漸近展開均質化法の適用により得られた微視的関係式 (3.18) 及び巨視的平衡式 (3.19) を有限要素法により近似表示する．
まず，前節にて導出した微視的関係式 (3.18) のマトリックス表記を以下に示す．

∫

(

T

T

∫

)

(

δ ε̇ Lχ,y + δq Qχ,y(q) dY =Y′)δ ε̇T L + δq T QR dY(3.20)Y

ここで L = D − F である．各マトリックスは次のように表せる．

(

)T

ε=

ε11 ε22 ε33 2ε12 2ε23 2ε31

χij
,y =

(

(
,

χ,y =

ij
ij
ij
ij
ij
χij
(11) χ(22) χ(33) χ(12) χ(23) χ(31)

)T
,

(
q=

χ11
χ22
χ33
χ12
χ23
χ31
,y
,y
,y
,y
,y
,y

χkl
(ij) =)T,

v1,1 v2,2 v3,3 v1,2 v1,3 v2,1 v2,3 v3,1 v3,2









Q=








σxx



1

 0

 0
R=
 0


 0
0
(
χ,y(q) =

0
σyy

0
1
0
0
0
0


σxy σxz 000000 σyx σyz 0
0 


0
0
0
0 σzx σzy 
σyy σyz 0
0
0
0 

σzz 0
0
0
0 
,
σxx σxz 0
0 


sym.
σzz 0
0 

σxx σxy 
σyy
0
0
σzz

T
0 0
0
0
0
0
0

0 0
0
0
0
0
0 

1 0
0
0
0
0
0 
 ,
0 1/2 0 1/2 0
0
0 


0 0
0
0 1/2 0 1/2 
0 0 1/2 0
0 1/2 0

22
33
12
23
31
χ11
,y(q) χ,y(q) χ,y(q) χ,y(q) χ,y(q) χ,y(q)

27

)
,

1
2

(

∂χkl
j
+
∂yj
∂yi

∂χkl
i

)
,

)
,

(

χij
,y(q) =

χij
(11)

χij
(22)

(
v=





F =





χij
(33)

χij
(12)

χij
(13)

χij
(21)

χij
(23)

χij
(31)

χij
(32)

)T
,

)T
v1 v2 v3

,

2σxx

0
0
2σzz

0
2σyy

σxy
σxy
0
(σxx + σyy ) /2


σxz


0


σzx
.

σzy /2
σxy /2(σxx + σzz ) /2

0
σzy
σzy
σzx /2(σzz + σyy ) /2

sym.

さらに，要素内の任意の点における変位速度 v 及び特性変位関数 χ をそれぞれ，要素の節点の変位速度 d 及び特性変位 χ(d) と形状関数 Ψ との線形結合によって，次のように表示する．形状関数 Ψ の具体形について，[付録 B] を参照されたい．

(
v = Ψḋ, ḋ =(
χ = Ψχ(d) , χ =(χ(d) =

χ11
(d)

χ22
(d)

T

T

ḋ1

χ33
(d)

χ12
(d)

T

χ23
(d)
(

χ31
(d)

(

T

, ḋN =

χ11 χ22 χ33 χ12 χ23 χ31χij(d)N =T)T

T· · · ḋN

ḋ2

)

, χij
(d) =

)

(

d˙N1 d˙N2 d˙N3

, χij =

T
χij
(d)1

ij
ij
χij
(d)N1 χ(d)N2 χ(d)N3

(

)
,

χij
χij
χij
1
2
3

T
χij
(d)2

···

)
,

T
χij
(d)N

)T

)
.

T

ここで，ḋN ，χij(d)N はそれぞれ，要素内 N 番節点の変位速度成分，特性変位成分である．また，要素内のひずみ速度 ε̇，変位速度勾配 q ，特性変位の偏微分 χ,y は節点変位速度 d 及び特性変位 χ を用いてそれぞれ次のように表すことができる．

ε̇ = B ḋ, q = E ḋ, χ,y = Bχ(d) , χ,y(q) = Eχ(d)(3.21)

ここで，B ，E は形状関数 Ψ を用いて表されるマトリックスであるが，その具体形については [付録 B] を参照されたい．式 (3.21) を式 (3.20) に代入することにより，微視構造における一つの要素に対する微視的方程式が得られ，次のようになる．

28

,

δ ḋ

T

[∫

(

∫

)

B LBχ(d) + E QEχ(d) dY −
T

T

(

T

T

)

]

B L + E QR dY=0(3.22)YY

このとき，任意の δ ḋ に対し式 (3.22) が成立するためには，次式が常に成立しなければ
ならない．

[∫

(

T

T

]

)

∫

( T
)
B L + E T QR dY

B LB + E QE dY χ(d) =Y(3.23)Y

つづいて，巨視的平衡式に移る．次式に巨視的平衡式 (3.19) のマトリックス表示式
を示す．

∫

(

T

H

T

H

)

∫
δv T Ṗ dS,

δε L ε + δq τ q dV =
Ω
∫
(
)
1
H
L − Lχ,y dV,
L =|Y | Y
∫
(
)
1
H
τ =
Q − Qχ,y(g) dV|Y | Y

St(3.24)

ここで，Ṗ ，χ,y(g) の具体形を以下に示す．

(
Ṗ =
(
χ,y(g) =)TṖ1 Ṗ2 Ṗ3

,

22
33
12
13
21
23
31
32
χ11
,y(q) χ,y(q) χ,y(q) χ,y(q) χ,y(q) χ,y(q) χ,y(q) χ,y(q) χ,y(q))

また，χ,y(g) は節点の特性変位 χ(d)g を用いて，次式で表せる．

(
χ,y(g) = Eχ(d)g , χ(d)g =

22
33
12
13
21
23
31
32
χ11
(d) χ(d) χ(d) χ(d) χ(d) χ(d) χ(d) χ(d) χ(d))次に，巨視的平衡式 (3.24) に式 (3.21) を代入することにより，全体構造における一つの要素に対する巨視的平衡式が得られ，次のようになる．

T

[∫

(

)

B L B + E τ E dV ḋ −
∫
(
)
1
H
L − LBχ(d) dV,
L =|Y | Y
∫
(
)
1
H
τ =
Q − QEχ(d)g dV|Y | Y
δ ḋ

T

H

T

H

Y

29

]

∫
T

Ψ Ṗ dS = 0,St

(3.25)

このとき，式 (3.25) が任意の δd に対して成立するには，次式が常に成立する必要がある．

K ḋ = f t ,
∫
( T H
)
K=
B L B + E T τ H E dV,∫Ωft =ΨT Ṗ dS

(3.26)

St

この式は要素の剛性方程式を表している．これを各要素について求め，全ての節点について重ね合わせると全体の構造剛性方程式を得ることができる．得られた構造剛性
方程式に境界条件を導入し，未知節点変位速度と未知節点力速度を決定する．それらからひずみ速度や応力速度などの各量が求められる．

30

第4章
多粒子モデルによるシリカ充填ゴムの力学特性評価
本章では，図 1.1 に示すようなシリカ粒子分散構造の力学応答を均質化法を用いて解析するが，周囲のゴム相及びシリカ粒子の界面に存在するゲル相は，それぞれの内部の絡み点の不均一分布によって力学特性が変化することが報告されている (6) ．そこ
で，まずゴム単相及びゲル単相それぞれの変形応答について，均質化法を用いて絡み点を不均一分散させた系の応答評価を行い，その結果を 2 章で示したモデル (図 2.3) にパラメータフィッティングすることで平均化して考慮する．上のスケールであるシリ
カ充填ゴム内部のランダムな粒子分散形態の均質化モデルでは，この平均化した構成式を導入する．

4.1

解析モデル

ゴム材に硫黄などの架橋剤を添加しゴム材を硫黄架橋した場合，図 4.1 に模式的に示すように，系内部に架橋点の多い（架橋点間のセグメント数 N が小さい）領域と少ない（セグメント数 N が大きい）領域が不均一に存在し，その大きさも場所によって大きく異なることが最近の実験的な研究により分かっている (6) ．このような硫黄架橋によるゴム内部の分子鎖の絡み点の不均一性に関して，その分布形態を観察する研究(7)が進められているが，実際のゴム内部における「セグメント数」の分布の詳細は明らかにされていない．そのため非晶性高分子材料での検討 (63) にならい，図 4.2 に示すようにゴム単相，ゲル単相での解析を，均質化法における微視領域において格子メッ

31

Fig.4.1 Schematic view of heterogeneous distribution of molecular chain(6) .

シュで区切り，各メッシュのセグメント数 N をランダムに変化させてその応答を均質化法により評価する．セグメント数の分布は図 4.3 に示すように (a) 正規分布，(b) ダブルピーク型分布とする．

4.1.1

ゴム相及びゲル相単相の絡み点不均一性に関する均質化法解析

図 4.3 において横軸は架橋点間の分子鎖セグメント数 N ，縦軸はユニットセル内でそのセグメント数を割り振った要素数 M である．ユニットセル内の平均初期セグメン
ト数 N0 = ΣN/ΣM はゴム相及びゲル相それぞれで定めるものとし，その N0 に基づいて分散範囲を 0.75N0 から 1.25N0 としたものを ND1 および DP1，0.50N0 から 1.50N0としたものを ND2 および DP2，0.25N0 から 1.75N0 としたものを ND3 および DP3 としている．このようなセグメント数 N の不均一分布に従い，ゴム相及びゲル相単相のユニットセルの要素に異なる初期セグメント数 N をランダムに割り当てたものを解析対象とする．
ただし，これまでの検討で最も不均一性の影響が大きくなったダブルピーク型分布

DP3 の応答をスケールアップモデルとして採用したので以降では DP3 の条件のみ記述
する．ゴムマトリクス相の平均初期セグメント数は N0 =12.0，ゲル相の平均初期セグ
メント数は N0 =8.0 として図 4.3 に示したようにユニットセル内にセグメント数をラン
ダムに割り振り，このユニットセルが周期的に繰り返されるゴム相及びゲル相単相の

32

Ν
21.0
17.4
13.8

L0

x2

10.2
6.6

y2

x1

0

0

3.0

y1

Unit cell(a) Rubber matrix

Ν
14.0
11.6
9.2

L0

x2

6.8
4.4

y2

0

x1

0

2.0

y1

Unit cell(b) Gel phase

Fig.4.2 Simulation model for non-uniform distribution of segment number N .

応答を，平面ひずみ状態，ひずみ速度 ϵ̇ = 10−1 ，最大引張り時のストレッチを λ2 = 4.0
として 1 サイクルの変形挙動を解析した．温度は変形過程を通して一定で T = 296[K]とし，各相の物理・化学架橋点の重み関数はどちらも w = 0.35 とした．その他の材料
定数は表 4.1 に示す通りである．各要素ごとに割り振られた初期セグメント数 N の値によって CαR ，CβR ，CBR の値は変化するため表 4.1 に示すのはユニットセル全体で平均した値である．また，これらの平均値をパラメータとする絡み点分布の均一なゴム相・ゲル相の解析 (均質化法を用いないものと同等) との比較を行う．

33

500

Number of elements M

Number of elements M

500

400

300

200

100

0

400

300

200

100

0

0.25N0 0.50N0 0.75N0 N0 1.25N0 1.50N0 1.75N0

0.25N0 0.50N0 0.75N0 N0 1.25N0 1.50N0 1.75N0

Segment Number N

Segment Number N(a) Normal distribution(b) Double peaks distribution

Fig.4.3 Distribution pattern of segment number ((a)normal distribution,(b)double peaks distribution).

34

Table 4.1 Average parameters in heterogeneous distribution of segment number
in Fig.4.2.
Rubber matrix
CαR (= nα kB T )
0.275
C1γ
3.0 × 102CβR = (nβ kB T )
0.40
C2γ
−0.50

CBR (= nB kB T )
0.50
mγ
6.4

Nα
12.0
C1D
3.0 × 103

Nβ
12.0
C2D
−0.50

NB
12.0
mD
4.8

Nα
8.0
C1D
2.5 × 103

Nβ
8.0
C2D
−0.50

NB
8.0
mD
4.8

Gel phase
CαR (= nα kB T )
0.405
C1γ
2.5 × 102

CBR (= nB kB T )
0.758
mγ
6.4

20

Homogeneous model
Heterogeneous model

True Stress Σt22[MPa]True Stress Σt22[MPa]20CβR = (nβ kB T )
0.599
C2γ
−0.50

15

15

10

10

5

0

Homogeneous model
Heterogeneous model

1

2

3

5

0

4

1

2

3

4

Stretch λ2(b) Gel phase

Stretch λ2(a) Rubber matrix

Fig.4.4 True stress-stretch relations (a)rubber matrix , (b)gel phase.

図 4.4(a) にゴム相，(b) にゲル相単相の解析結果の真応力-ストレッチ関係を，絡み点の不均一性を考慮した場合と均一のそれを比較して示す．ゴム相，ゲル相ともに絡み点に不均一分布を導入すると負荷変形開始直後から応力が上昇している．これはユニットセル内において割り振られたセグメント数が小さい要素 (硬い要素) では変形が起こりにくく，それを補うようにセグメント数が大きい要素 (軟らかい要素) に変形が35

集中する，という現象がゴムの内部で発現することに起因する．この現象はシリカ粒子等のゴムよりも剛性の非常に大きいものを充填した時と同じ現象であり，変形が集中する剛性の小さい部分において応力が非常に高くなり，配向硬化を引き起こすというものである．この効果が絡み点不均一分布の導入により得られたため，図 4.4 に示すように応力が上昇したと考えられる．

Table 4.2 Parameters for constitutive equations Fig4.5 (a)rubber matrix, (b)gel
phase ﬁtted to the homogenized FEM results of Fig.4.2.
Rubber matrix
CαR (= nα kB T )
0.324
C1γ
6.0 × 102CβR = (nβ kB T )
0.479
C2γ
−0.50

CBR (= nB kB T )
0.606
mγ
6.4

Nα
10.0
C1D
6.0 × 103

Nβ
10.0
C2D
−0.50

NB
10.0
mD
4.8

Nα
7.0
C1D
8.0 × 103

Nβ
7.0
C2D
−0.50

NB
7.0
mD
4.8

Gel phase
CαR (= nα kB T )
0.463
C1γ
8.0 × 102CβR = (nβ kB T )
0.685
C2γ
−0.50

CBR (= nB kB T )
0.886
mγ
6.4

この不均一分布による効果を均一なゴム相・ゲル相の構成式に加味するため，ゴム相，ゲル相ともにひずみ速度を ϵ̇ = 10−1 ，ϵ̇ = 10−5 ，最大引張り時のストレッチを λ2 =2.0，

3.0，4.0 と変化させた時の真応力-ストレッチ関係を均質化法により求め，パラメータ
フィッティングを行った．図 4.5 にゴム相，図 4.6 にゲル相の各ひずみ速度における真応力-ストレッチ関係をそれぞれ示す．フィッティングにより得られたパラメータを表4.2 に示す．

36

Heterogeneous model(DP3)Pseudo-heterogeneous model

8

6

4

2

0

1

2

3

Heterogeneous model(DP3)Pseudo-heterogeneous model

10

True Stress σt22[MPa]True Stress σt22[MPa]

10

8

6

4

2

0

4

1

2

3

Stretch λStretch λ(a) Strain rate = 0.00001(b) Strain rate = 0.1

4

Fig.4.5 True stress-stretch relations of rubber matrix by simulation.

Heterogeneous model(DP3)Pseudo-heterogeneous model

15

15

10

10

5

0

20

Heterogeneous model(DP3)Pseudo-heterogeneous model

True Stress σt22[MPa]True Stress σt22[MPa]

20

1

2

3

5

0

4

1

2

3

Stretch λStretch λ(a) Strain rate = 0.00001(b) Strain rate = 0.1

4

Fig.4.6 True stress-stretch relations of gel phase by simulation.

37

4.1.2シリカ粒子分散構造の均質化法モデル

Silica particle.

Rubber matrix.

Gel phase.

50nm(a)Tire.

(b)Silica filled rubber.

(c)Bunching structure.

Fig.4.7 Structure of silica ﬁlled-rubber(4) .

(a) All silica particles(b)Proximal particles

Fig.4.8 Distribution of silica particles Red colored particles are proximal particles with distance of smaller than 20nm(5) .

シリカ充填ゴムは図 4.7(b) に示すようにシリカ粒子をランダムに分散した状態で含有しており，さらに拡大すると (c) に示すようにシリカ粒子がゲル相を介して数珠繋ぎ状に連結していることが観察されている (4) ．また図 4.8(a) はシリカ充填ゴム内部のシリカ粒子にのみ色を付けて 3 次元的に可視化したものであり，粒子間距離が 20nm 以

38

下の接近粒子を赤，それ以外の粒子を青で色付けしている (5) ．そして赤で示す接近粒子のみを抽出したものを図 4.8(b) に示す．この図からもシリカ粒子がランダムに分散している状態が確認出来，また粒子間距離が非常に近いもの，つまり図 4.7(c) のようにゲル相を介して粒子が連結しているような部分も多数存在していることが示唆されている．そこで本章では平面ひずみ条件下においてランダムな粒子分散形態で，かつゲル相を介して数珠繋ぎ状に粒子が連結する構造を有する均質化有限要素解析の対象

x2

L0

とする．

x1

0

y2
0

y1

Unit cell

x2

0

L0

(a) Grid points model

x1

y2
0

y1

Unit cell(b) Random modelFig.4.9 Simulation model(silica ﬁlled-rubber).

図 4.9 に解析モデルの例を示す．青で着色した部分がシリカ粒子，粒子周りの黄緑がゲル相，それ以外の赤がゴムマトリクス相である．シリカ，ゴム，ゲル相はボクセ
ル要素毎に決めており，平面ひずみ条件下で y1 ，y2 方向にそれぞれ 100 分割のボクセ
ル分割を行っている．シリカ粒子の直径は全て同じとし，検討する物性としてシリカ

39

粒子の数と含有率，ゲル相の厚さ等を変えた解析を行う．シリカ粒子中心の配置はユニットセル内のボクセル分割により出来た格子点上 (ユニットセル境界部分含む) からランダムに選択して配置するが，(a) シリカ粒子の並びが y1 ，y2 方向に常に平行となるように，ボクセルメッシュより大きいグリッド上に配置した規則格子モデル，(b) 規
則格子モデルからさらに乱数を用いて粒子を偏心させたランダムモデル，の 2 種類を考慮した．粒子数は 36 個，粒子含有率は 20%，ゲル相厚さはユニットセル内のゲル相割合が 8.8%となるように設定した．
解析には，周期的微視構造を有する材料全体を表現する座標系 xi と微視構造を表現する座標系 yi の二変数を用い，変位を漸近展開する均質化理論 (62) に基づき定式化した有限要素法 (59) を用いる．本研究で直接用いる均質化理論，並びに，有限要素方程式の具体形，計算手順の詳細については 3 章及び文献 [5,15] を参照されたい．
巨視的に一様な単軸変形を与えるものとして，巨視座標系の x2 方向に u̇ = 100[mm/min]の一定な変形速度で，最大伸びが λ2 = 1.5 になるまで引張変形を与えた後除荷を行い，1 サイクルの変形挙動を解析する．ゴムの非圧縮性を満足させるためのペナルティ定
数は ṗ = 50 とする．シリカ粒子の剛性は，ゴム材の剛性に比べて十分大きいと考え，
ゴムマトリクス部とゲル相の計算の安定性と結果にほとんど影響を与えない値として，
縦弾性係数 E = 100[MPa], ポアソン比 ν = 0.3 とした．材料の温度はゴム相，ゲル相単相の解析と同様，変形過程を通して一定で T = 296[K] とした．
本モデルについて，次節にてモデルの評価精度についてユニットセルの等方性や解像度の角度から議論し，次々節にて本モデルを用いたシリカ充填ゴムの力学特性評価
を行う．

40

4.2

解析モデルの検証

本節では 4.1.2 で示したシリカ充填ゴムの解析モデルについて，ユニットセルの等方性，解像度の影響を調べ，シリカ充填ゴムの力学特性評価の精度について検討する．

4.2.1

モデルの等方性

実際のシリカ充填ゴムはランダムな粒子配置を呈しているが，引張り方向による変形応答の違いというものはなく，概ね等方性が満たされている材料と考えられる．そこで本モデルにおいても実際のゴムと同様に等方性が満たされているか否かを図 4.9の各モデルにおいて検証する．

Vertical Tension
Horizonal Tension

True Stress Σt22[MPa]

4

3

2

1

0

1

1.1

1.2

1.3

1.4

1.5

Stretch λ2(b)True Stress-Stretch relations(a)Unit cell

Fig.4.10 (a)Unit cell model and (b)true stress-stretch relations of vertical tension and horizonal tension by grid points model.

始めに規則格子モデル (図 4.10(a)) に矢印で示す方向に引張りを与えた場合の真応力ストレッチ関係を図 4.10(b) に示す．両者はほとんど一致していることから，図 4.9(a)の規則格子モデルの応答は等方性を満たしている．
同様にランダムモデル (図 4.11(a)) に矢印で示す方向に引張りを与えた場合の真応

41

Vertical Tension
Horizonal Tension

True Stress Σt22[MPa]

4

3

2

1

0

1

1.1

1.2

1.3

1.4

1.5

Stretch λ2(b)True Stress-Stretch relations(a)Unit cell

Fig.4.11 (a)Unit cell model and (b)true stress-stretch relations of vertical tension and horizonal tension by random model.

力-ストレッチ関係を図 4.11(b) に示す．規則格子モデルと比べるとわずかに方向による差を生じているが，ほとんど一致していることから，図 4.9(b) のランダムモデルも概ね等方性を満たしている．

4.2.2

解像度の検証有限要素法は，本来連続であり細分化は無限 (=原子 1 個 1 個に至るまで) である物
質を，有限の要素で近似的に代表させることで物質に生じる変位や応力などを数値計
算によって評価する方法であり，要素数が多いほど，より精度の高い評価が出来るものとされている．そこで本節では図 4.10(a) の規則格子モデルについて，要素分割数が

y1 ，y2 方向にそれぞれ 100 分割，つまり crossed triangles 要素で 10000 要素のものと，粒子の配置，直径，含有率は同じで分割数のみ 2 倍の 200 分割，crossed triangles 要
素で 40000 要素のモデルによる解析結果を比較することによりモデルの解像度を検討する．
図 4.12(a) に分割数 100，(b) に分割数 200 のユニットセルを，(c) に真応力-ストレッチ関係を示す．全体的な応答は概ね一致しているものの，分割数 200 モデルによる結42

meshing 100
meshing 200

True Stress Σt22[MPa]

4

3

2

1

0

(a)meshing 100 × 100(b)meshing 200 × 200

1

1.1

1.2

1.3

1.4

Fig.4.12 Comparison of true stress-stretch relations by (a) meshing 100 × 100
and (b) meshing 200 × 200.

果の方が負荷時，除荷時ともに高い応力を示しており，1 サイクルのヒステリシスルー
プが上にシフトしている．この要因に関して最大引張り時 (λ2 =1.5) における分子鎖ス
トレッチ λc の分布及び引張り方向の応力 σ22 の分布から考察する．
図 4.13 に最大引張り時における分子鎖ストレッチ λc の分布を，図 4.14 に引張り方向の応力 σ22 の分布をそれぞれ示す．λc の分布について，(a)(b) どちらのモデルにおいても粒子が密集している部分に変形が集中している傾向が見られ，粒子の不均一分散による内部応力状態には劇的な変化はない．しかしながら，ユニットセル中央の粒子がゲル相を介して連結している部分や，ユニットセル右端の粒子付近に着目すると，
分割数 200 のモデルの方が局所的な変形をよくとらえており，それにより λc が大きくなっていることがわかる．それに対応して σ22 が著しく大きくなっている部分がある
ことがわかる．このことから分割数の大きいモデルの方が変形集中部の λc 及び σ22 を表現し，それが図 4.12(c) の真応力-ストレッチ関係に現れたものと考える．
図 4.15 に有限要素メッシュを各相ごとに色分けしたものを示し，特にユニットセル
内の黒枠の部分を拡大して右に示す．拡大図の粒子連結部のゲル相に着目すると，(a)の分割数 100 の場合はシリカ粒子を接続しているゲル相のメッシュが不足していることがわかる．
以上の結果から，巨視的な (平均化された) 系の応答については分割数による違いは小さいが，内部の粒子連結部における変形集中の精度はメッシュ数に依存することが

43

1.5

Stretch λ2(c)True Stress-Stretch relations

わかった．しかしながら，200×200 の計算は計算機コストや数値計算上の安定性などの課題があるため，以降ではシリカ充填ゴムの巨視的な力学特性は分割数の粗いモデ
ルにおいて評価し，粒子周りの変形集中部の変形の詳細は，対象とする領域を小規模なセルでモデル化し直して検討することとする．

44

λc
2.44

2.15

1.86

1.58
1.29
1.00
(a)meshing 100 × 100(b)meshing 200 × 200

Fig.4.13 Distribution of chain stretch λc at the maximum system stretch of λ2 = 1.5

σ22[MPa]
18.39

13.74

9.10

4.46
-0.19
-4.83

(a)meshing 100 × 100(b)meshing 200 × 200

Fig.4.14 Distribution of tensile stress σ22 in the loading direction at the maximum system stretch of λ2 = 1.5

45

(a)meshing 100 × 100(b)meshing 200 × 200

Fig.4.15 Deformed ﬁnite element mesh at the maximum system stretch of λ2 = 1.5

46

4.3

シリカ充填ゴムの力学特性評価

前節における解析モデルについての検証結果から，分割数 100 で，等方性が満たされている規則格子モデルとランダムモデルのうち，より現実の粒子配置に近いランダ
ムモデルを用いてシリカ充填ゴムの巨視的応答の評価を行う．本節では，最大引張り
時 (λ2 =1.5) における応答，また 1 サイクルの変形過程における応答を検討することによりシリカ充填ゴムの力学特性を評価する．

4.3.1

最大引張り時における応答

図 4.16(a) に真応力-ストレッチ関係を，最大引張り時について，(b) 有限要素メッシュ
図を各相ごとに色分けしたもの，(c) 分子鎖ストレッチ λc の分布，(d) 引張り方向応力

σ22 の分布を合わせて示す．(a) の真応力-ストレッチ関係より負荷から除荷に移行しても急激な応力の低下は見られず，ヒステリシスロスも小さい．

(c) の分子鎖ストレッチ λc 分布について，図 4.16(b) の (1) で示した部分のように周囲を剛性の高い粒子に囲まれた狭い領域では，ユニットセルの延伸に伴い非常に大きな変形が発生している．また図 4.16(b) の (2) で示したように引張り方向に並んだ粒子と粒子の間の部分でも，粒子と粒子を繋ぐように λc の値が大きいラインが確認出来る．
逆に図 4.16(b) の (3) のように粒子が凝集しているところでは，構造の「凹み」部分はほとんど変形せず，その外側のゴム相に大きな変形が見られることから粒子のクラス
ターが 1 つの介在物のようにふるまっている．
次に (d) の引張り方向応力 σ22 分布に着目すると，図において赤や黄色で示されている応力値が高い部分は，粒子を引張り方向に連結するゲル相部分である．この部分は
図 4.16(c) の λc の分布では特筆するほど大きい値を示してはいない．これはゴム相に比べゲル相の配向硬化が，小さいストレッチにおいて発生することに起因する．また黄緑色の比較的大きな応力がセル上下を支持するように分布しているが (図 4.16(b) の(2) や (3) 付近)，これは λc が比較的大きいゴム相に対応している．さらに，青色はカラーバーに示したように引張り方向の応力が負となった部分の存在を示している．これは図 4.16(b) の (4) に示したように複数の粒子に囲まれたゴム相部分である．この部
分では，引張りによって左側の粒子群が右側のそれに対して相対的に下に変位するが，この際に粒子間で σ22 方向の圧縮を生じたものと考えられる．

47

True Stress Σt22[MPa]

4

3

(4)
(1)

2

(3)
1

(2)
0

1

1.1

1.2

1.3

1.4

1.5

Stretch λ2
(b)

(a)

σ22[MPa]

λc
2.44

18.39

2.15

13.74

1.86

9.10

1.58

4.46

1.29

-0.19

1.00

-4.83
(d)

(c)

Fig.4.16 (a)True stress-stretch relations and deformation behaviors under 1cycle loading with maximum stretch of λ2 =1.5,(b)Deformed ﬁnite element mesh, (c)distributions of chain stretch, (d)distribution of tensile stress.(random model)

4.3.2

1 サイクルの変形過程における応答

図 4.17 に負荷 → 除荷過程における分子鎖ストレッチ λc の分布の変化を，上側 (a)に負荷時，下側 (b) に同じストレッチの除荷時におけるものを並べて示す．なお，本48

(a) Loading

λc
2.44
2.15
1.86
1.58

λ2=1.1

1.29(b) Unloading

λ2=1.2λ2=1.3λ2=1.4λ2=1.5

1.00

Fig.4.17 Change in distribution of chain stretch.(random model)

シミュレーションでは 1 サイクルの変形過程において最大引張り状態での保持時間を0 としているため，λ2 =1.5 における負荷時と除荷時の分布 (右端の上下の図) は全く同じ分布であることに注意されたい．一見すると負荷時と除荷時の λc の分布に差異は見られないが，λ2 =1.2 や λ2 =1.3 の図からわかるように，粒子の集合体付近のゴム相など局所的に変形が集中していた部分においては，若干ではあるが除荷時の方が λc が大
きい (コントラストが明確) であることがわかる．
次に図 4.18 に引張り方向応力 σ22 の分布の変化を，上側 (a) に負荷時，下側 (b) に除荷時のものを並べて示す．こちらも負荷と除荷で全く同じに見えるが，例えば λ2 =1.2の図で比較すると，除荷時の方がコントラストがぼんやりしている．これは先の応力ストレッチ関係で系の平均応力が負荷時より除荷時の方がわずかに小さいことに対応

49

(a) Loading

σ22[MPa]
18.39
13.74
9.10
4.46

λ2=1.1λ2=1.2λ2=1.3λ2=1.4λ2=1.5(b) Unloading

-0.19
-4.83

Fig.4.18 Change in distribution of tensile stress.(random model)

する．これは負荷変形時に発生する非アフィン変形により 1 分子鎖あたりのセグメン
ト数 N が増加し，変形抵抗が減少したためである．これにより図 4.16(a) の真応力-ス
トレッチ関係においてヒステリシスロスが発生している．
以上の検討から，負荷過程において変形集中が見られるのは粒子周りのゴム相やゲル相であり，その顕著な応力上昇が系の巨視的な応答を支配していることがわかる．次章ではこの変形集中について小規模モデルによる詳細な検討を行う．

50

第5章
シリカ充填ゴムにおける変形集中機
構の検討
4 章では多数のシリカ粒子をランダムに配置した，実際の内部構造を模擬した解析モデルによりシリカ充填ゴムの巨視的な応答の評価を行ったが，個々の粒子近傍での力学状態ならびに境界条件がシリカ粒子周りの変形集中で複雑となり，シリカ-ゲル相の変形集中の本質に迫ることが難しい．議論するためには，変形集中部の解像度を十分に上げる，かつ問題の切り分けのためにシンプルかつ極端な境界条件とした方が新
たな知見を得られるものと考える．そこで本章では，ユニットセル内に 2 粒子のみ含有させた小規模なモデルにおいて，シリカ粒子の数珠繋ぎ構造による応答変化について議論する．

5.1

解析モデル

シリカ充填ゴム内部の粒子の数珠繋ぎ構造を構成する最小単位として，平面ひずみ
条件下で 2 粒子が引張り方向に直鎖状に繋がったモデル (図 5.1(a))，及び，引張り方
向に対して 45◦ 交互に傾いて繋がったモデル (図 5.1(b)) を用いてシリカ充填ゴムの粒
子周りの局所的な変形応答の評価を行う．

4 章とは異なりユニットセル内のシリカ粒子含有率は 30%，ゲル相厚さはユニット
セル内のゲル相割合が 12.54%としているが，粒子含有率に対するゲル相割合の比率は4 章のモデルとほとんど等しく，粒子径の大きさに対するゲル相厚さは 4 章と同程度のモデルである．また導入するゴム相及びゲル相の物性値は絡み点不均一性を考慮した

51

L0

x2

x1

0

y2
0

y1

x2

Unit cell

L0(a) 0° model

y2

x1

0

0

y1

Unit cell(b) 45° model

Fig.5.1 Simple simulation models to investigate the eﬀects of the direction ofsilica bunching structure in rubber.

表 4.2 に示すものを用い，その他解析条件は 4 章と同じく非圧縮性を与えるペナルティ
定数を ṗ = 50，シリカ粒子の縦弾性係数 E = 100[MPa]，ポアソン比 ν = 0.3，材料温度を T = 296[K] で一定とする．巨視的に一様な単軸変形を与えるものとして，巨視座
標系の x2 方向に u̇ = 100[mm/min] の一定な変形速度で，最大伸びが λ2 = 1.5 になるまで引張変形を与えた後除荷を行い，1 サイクルの変形挙動を解析する．

5.2

解析結果

図 5.2(a) に 0◦ モデルの解析結果の真応力-ストレッチ関係を示す．また最大引張り時の λ2 =1.5 における (b) 初期セグメント数からの変化量 ∆N の分布，(c) 分子鎖スト
レッチ λc の分布，(d) 引張り方向応力 σ22 の分布，(e) 回転 θ の分布を示す．なお (d)の σ22 について，粒子連結部のゲル相の一部のみが σ22 =86.0[MPa] 程度まで上昇しており，その他の相は大きいところでも 15.0[MPa] 程度である．この要因については後52

述するが，これにより σ22 の最大値にスケールを合わせると σ22 の詳細な分布状態を確認出来ないためスケールの最大値を σ22 =16.0[MPa] とし，この値を超えた部分は全て赤で示している．4 章で用いた解析モデルとは粒子含有率やゲル相厚さが異なること，

0◦ モデルや 45◦ モデルは粒子配置が特別なユニットセルであることから 4 章の結果と値の大きさの比較は行わず，その特徴を検討する．そのため (a)(c)(d) のスケールは 4 章のものとは異なることに注意されたい．

0◦ モデルによる解析では負荷過程において応力が軟化した後，更に引張りを与える
ことで配向硬化が見られ，4 章の巨視的な応答では見られなった 3 次関数的な曲線を描いている．また負荷反転後に応力の大幅な低下が見られ，それに伴いヒステリシス
ロスも大きく発生していることがわかる．
次に図 5.2(c) の λc 分布，(d) の σ22 分布，(e) の θ 分布から 0◦ モデルの詳細な変形応答について言及する．0◦ モデルはシリカ粒子が引張り方向に直線的に連結されており，剛性の高いシリカ粒子はほとんど変形しないためユニットセルの中心軸付近では系全体に与えた引張りが粒子連結部のゲル相に集中している．加えて，(e) の回転 θ の分布より，引張りに伴い粒子周りのゲル相はそれぞれ粒子側面方向に回転しているため，より大きな変形が中心軸付近のゲル相に加わっている．そのため引張り方向応力 σ22 が16.0[MPa] 以上の高い値を示している．さらに，Poisson 収縮によって 2 つの粒子の側面部分の間隔が縮小されるので，その間にあるゴム相は，他のゴム相に比べ大きな λc ，

σ22 を生じている (図 (c),(d) の横方向の帯)．
次に図 5.3 に負荷 → 除荷過程における分子鎖ストレッチ λc ，図 5.4 に引張り方向応力 σ22 それぞれの分布の推移を，上側 (a) に負荷時，下側 (b) に除荷時のものを並べて示す．4 章と同じく右端の λ2 =1.5 における負荷時と除荷時の分布は同じである．負荷過程においてまず粒子周りのゲル相の λc が周囲に比べて大きくなり，λ2 =1.3 近傍から粒子連結部のゲル相の変形が顕著になっていく様子，それに伴う σ22 の顕著な上昇が確認出来る．λc の最大値は 2.80 と非常に大きく，これは引張り方向の主ストレッチ

λ2 に換算すると 4.75 程度に相当する．4.1 節の図 4.6 で示した，解析モデルに導入しているゲル相の応答から，負荷過程 λ2 =3.0 近傍から配向硬化が見られることがわかる．この 0◦ モデルにおいても λ2 =3.0 を大幅に超える変形がゲル相に加わっているため，図 5.2(a) の真応力-ストレッチ関係においてもこのゲル相の影響により配向硬化が
見られた．

53

一方除荷については負荷過程と全く同じ順序で分子鎖ストレッチ λc が元の値に収束しており，引張り方向応力 σ22 については 4 章の巨視的な応答と同様負荷時よりも応力値が小さいが，特に粒子連結部のゲル相の応力値が大幅に小さくなっていることがわ
かる．図 5.2(b) のセグメント数変化量 ∆N の分布から，粒子連結部のゲル相の ∆N は

2.9 程度，他のゲル相は 1.5 程度と λc の大きさに対応してセグメント数が増加しており，それにより変形抵抗が減少するため，応力の低下が顕著になったと考えられる．このことからシリカ充填ゴムの 0◦ モデルに発生するヒステリシスロスは，特に非アフィン変形の顕著な粒子連結部において大幅に応力が低下することに起因することがわかった．

54

ΔN
4.84
6

3.87
True Stress Σt22[MPa]

5

2.91

4
3

1.94
2

0.97

1
0

λc

1

1.1

1.2

1.3

1.4

0.00

1.5

Stretch λ2
(a)

(b)

σ22[MPa]

θ[°]

2.80

16.00

2.44

12.66

2.08

9.32

1.72

5.98

1.36

2.64

1.00

-0.70

62.8

37.7

12.7

-12.4
-37.5
-62.5
(d)

(c)

Fig.5.2 (a) True stress-stretch relations and deformation behaviors for an
average stretch λ2 =1.5 distribution of (b)changes of average numberof segments, (c) chain stretch and (d) tensile stress. (0◦ bunching
structure model)

55

(e)

λc (a) Loading
2.80

2.44

2.08

1.72

λ2=1.1(b) Unloading

λ2=1.2λ2=1.3λ2=1.4λ2=1.5

1.36
1.00

Fig.5.3 Change in distribution of chain stretch. (0◦ bunching structure model)

σ22[MPa] (a) Loading
16.00

12.66

9.32

λ2=1.1(b) Unloading

λ2=1.2λ2=1.3λ2=1.4λ2=1.5

5.98
2.64
-0.70
Fig.5.4 Change in distribution of tensile stress. (0◦ bunching structure model)

56図 5.5(a) に 45◦ モデルの結果をまとめて示す．比較のため各スケールは 0◦ モデルと合わせている．なお，45◦ モデルの場合は σ22 の値の最大値が 16.0[MPa] 程度である．

(a) の真応力-ストレッチ関係について，0◦ モデルとは異なり負荷過程における顕著な応力軟化や配向硬化は見られなかったが，除荷時の応力の低下によるヒステリシス
ロスの発現は確認出来る．図 5.5(c)，(d) 及び (e) から，45◦ モデルにおいても 0◦ モデルと同様にシリカ粒子間の距離が最小となる中心軸付近のゲル相に変形が集中し，顕著な応力上昇が生じていることがわかる．加えて，引張りに伴い粒子位置が中心軸方向へずれるため，(e) の回転分布で示されるように粒子の両側面で回転角度の大きさが異なる．このため粒子の両側面のゴム相の分子鎖ストレッチ λc は非対称となり，回転角度の大きい側面のゴム相の変形が顕著になる．これによりゲル相の変形，引張り方向
応力 σ22 の値も緩和されている．
最後に図 5.6 及び図 5.7 に負荷 → 除荷過程における λc ，σ22 それぞれの分布の変化を示す．0◦ モデルとは異なり，負荷過程では粒子連結部のゲル相だけでなく回転角度の大きい側面のゴム相が，同程度の λc を示しながら系全体が延伸している．このため

σ22 の大きい部分が粒子周りで S 字を描くように分布し，つり合いの条件から中心軸上の粒子にも大きな応力が分布している．最大引張り時における λc の最大値は粒子連結部のゲル相や回転角度の大きい側面のゴム相で 2.3 程度であるが，そこまで伸びている要素の範囲は非常に小さく，ゲル相では大きく延伸している側で 2.0 程度である．

λc =2.0 は引張り方向の主ストレッチ λ2 に換算すると 3.3 程度に相当する．4.1 節の図
4.5 及び図 4.6 から λ2 =3.3 程度ではどちらも配向硬化は生じない．そのため，図 5.5(a)の真応力-ストレッチ関係において 0◦ モデルのような顕著な配向硬化は見られなかったと考えられる．
一方，除荷については分子鎖ストレッチ λc ，引張り方向応力 σ22 ともに負荷時と形態は変わらず，コントラストが弱っただけとなる．ここで図 5.2(b) の 0◦ モデルのセグ
メント数変化量 ∆N と図 5.5(b) の 45◦ モデルの ∆N のカラーバーの範囲を比べると，

45◦ モデルの方が全体的にセグメント数増加量が大きいことがわかる．これは 0◦ モデル
ではゲル相への変形集中が顕著であったのに対して 45◦ モデルは粒子位置のずれによりゲル相に集中する変形をユニットセル全体に緩和したためである．これにより，45◦
モデルでは除荷時における粒子連結部のゲル相の応力低下だけでなく，ユニットセル
全体での応力低下によりヒステリシスロスが発現することが確認出来た．

57

ΔN
4.84

6

True Stress Σt22[MPa]

5

3.87
4

2.91

3
2

1.94
1
0

0.97
1

1.1

1.2

1.3

1.4

1.5

Stretch λ2

0.00

(a)
σ22[MPa]

λc

(b)
θ[°]

2.80

16.00

62.8

2.44

12.66

37.7

2.08

9.32

12.7

1.72

5.98

-12.4

1.36

2.64

-37.5

1.00

-0.70

-62.5

(d)

(c)

(e)

Fig.5.5 (a) True stress-stretch relations and deformation behaviors for an average stretch λ2 =1.5 distribution of (b)changes of average number ofsegments, (c) chain stretch, (d) tensile stress and (e) material rotation. (45◦ bunching structure model)

58

λc

(a) Loading

2.80

2.44

2.08
λ2=1.11.72

λ2=1.2λ2=1.3λ2=1.4λ2=1.5(b) Unloading

1.36
1.00
Fig.5.6 Change in distribution of chain stretch. (45◦ bunching structure model)

σ22[MPa] (a) Loading
16.00

12.66

9.32
λ2=1.15.98

λ2=1.2λ2=1.3λ2=1.4(b) Unloading

2.64
-0.70
Fig.5.7 Change in distribution of tensile stress. (45◦ bunching structure model)59

λ2=1.5

第6章
結論
シリカ充填ゴムは，シリカ粒子の分散形態や，粒子界面の結合状態を調整する結合剤の量を変えることで，これまでの主流である CB 充填ゴムより広範囲に力学特性を制御でき，さらには環境負荷が少ない．本研究ではシリカ充填ゴムの微視的な数珠繋ぎ構造をモデル化し，その力学応答を明らかにするため，非アフィン分子鎖網目理論
によるゴム・ゲル相の物理・化学架橋を考慮した構成式の構築，均質化法に基づく有限要素シミュレーションによる多数のシリカ粒子分散構造の力学応答評価や少数の数
珠繋ぎ構造での基礎物理メカニズムの解明などを行った．得られた結果を以下にまとめて示す．
第 2 章では，物理架橋点の解消を許容する非アフィン分子鎖網目理論により一般化
した粘弾性 8 鎖モデルに，化学架橋点の解消も考慮する複合型非アフィン分子鎖網目モデルへの拡張を行った．さらに，負荷から除荷へ反転する際に観察される大きな応力低下現象を表現する構成式を提案し，シリカ未充填のゴムの新たな構成式を提案した．
第 3 章では，2 変数漸近展開理論に基づく均質化法の基本的な考え方を述べた後，2章で示したゴム材の構成式を更新ラグランジュ法に基づく有限要素均質化法として定式化した．
第 4 章では，数珠繋ぎ構造を呈するシリカ充填ゴム内部のランダムな粒子分散形態を模擬したユニットセルモデルを構築し，ユニットセルの等方性，解像度の影響を検証した後，シリカ充填ゴムの巨視的な変形応答を評価した．シリカ粒子を含有する系の解析の前に，まず実験により報告されている絡み点の不均一性について，ゴム相及びゲル相それぞれの単相について均質化法で格子上に分割した要素ごとに，絡み点数を決定するセグメント数 N をランダムに変えた解析を行い，均一とした場合との差を60

構成式のパラメータで再フィッティングし直すことで考慮した．そのゴム相・ゲル相の構成式を用いて，36 個のシリカ粒子を配置させた系の均質化法有限要素シミュレー
ションを行い，シリカ粒子の並びが縦横整列している規則格子モデルとランダムモデ
ルの等方性を確認した後，現実の粒子配置に近い後者でシリカ充填ゴムの巨視的応答に及ぼす内部構造の変化について議論した．λ2 =1.5 までの繰り返し変形において，ユニットセルの延伸に伴い粒子に囲まれた狭い領域や引張り方向に並んだ粒子間のゴム相に大きな変形が集中していることや，粒子が凝集しているところでは粒子のクラス
ターが 1 つの介在物のようにふるまうことなどを明らかにした．さらに負荷変形時に生じる非アフィン変形，すなわち絡み点数の減少によって変形抵抗が減少し，ヒステ
リシスロスをもたらすことを示した．
第 5 章では，分散粒子間で生じる変形集中及び粒子周りの変形メカニズムを，より明確にするために簡素化した境界条件でユニットセル中に 2 粒子が引張り方向に直鎖状に繋がった 0◦ ，ならびに引張り方向に対して 45◦ に傾いて繋がった 45◦ モデルの 2 種類の極端な条件で引張りシミュレーションを行った．0◦ モデルの場合，粒子が引張り方
向に直線的に連結されているため粒子連結部のゲル相に変形が集中し，負荷変形後期における配向硬化や顕著な非アフィン変形を生じてヒステリシスロスが増大した．一
方 45◦ モデルは引張りに伴い粒子位置が中心軸方向へ相対的にずれるためゲル相への変形集中は緩和され，0◦ モデルに比べユニットセル全体での非アフィン変形を生じてヒステリシスロスを担っていることが示された．
以上のように，ゴム・ゲル相の物理・化学架橋点の解消や，負荷反転時の応力低下を考慮した構成式モデルの提案，絡み点の不均一性によるゴム・ゲル相の応答変化，多数のシリカ粒子分散構造におけるヒステリシス発現と局所変形メカニズムの議論，2粒子数珠繋ぎ構造でのヒステリシス発現メカニズムの検討などを行った．これにより，シリカ充填ゴムの変形に及ぼす微視的なメカニズムの理解に関わる基本的な情報が得られたと考える．シリカ充填ゴムを用いた製品の実設計においては 3 次元的な変形挙動を呈するので，本論文で解析対象とした平面ひずみ状態でのシミュレーション結果に比べ，変形の自由度は高くなり，より変形抵抗の小さい応答となることが考えられる．したがって，今後の材料の実設計では CB 充填ゴムの場合と同様の本研究成果を3 次元に拡張した大規模解析による力学特性の評価が望まれる．

61

参考文献
(1) 深堀美英, 設計のための高分子の力学, (2000), 技報堂出版.
(2) Mullins, L., Thixotropic Behavior of Carbon Black in Rubber, Rubber Chemistry
and Technology, Vol.23, (1948), pp.281-300.
(3) Bergstöm, J.S. and Boyce, M.C., Constitutive Modeling of the Large Strain TimeDependent Beahvior of Elastomers, Journal of Mechanical Physics and Solids,
Vol.46, No.5(1998), pp.931-954.
(4) 住友ゴム工業株式会社, TEM 画像.
(5) Courtesy of SRI R&D Ltd.
(6) 柴山充弘, 日本ゴム協会誌, 84-1,(2011),14-20.
(7) 中島健, 岩蕗仁, 浦部匡史, 伊藤万喜子, 藤波想, 西敏夫, 日本ゴム協会誌, 857,(2012),46.
(8) D. Wang, S. Fujinami, K. Nakajima, S. Inukai, H. Ueki, A. Magario, T.
Noguchi, M. Endo, T. Nishi, Visualization of nanomechanical mapping on
polymer nanocomposites by AFM force measurement, Polymer, Vol.51, (2010),
pp.2455-2459.
(9) 池田裕子, シンクロトロン X 線による架橋天然ゴムの伸長結晶化解析, 日本ゴム協会誌, 84-1,(2011),pp.29-36(10) Gao, J. and Weiner, J. H., Computer Simulation of Viscoelasticity in Polymer
Melts, Macromolecules, 25, (1992), 1348-1356.
(11) Ananyo, Bandyopadhyay, Pavan, K. Valavala, Thomas, C. Clancy, Kristopher, E.
Wise. and Gregory M. Odegard., Molecular modeling of crosslinked epoxy polymers: The eﬀect of crosslink density on thermomechanical properties, Polymer,
52, (2011), 2445-2452.
62

(12) Ying, Li, Martin, Kroger. and Wing, Kam, Liu., Primitive chain network study
on uncrosslinked and crosslinked cis-polyisoprene polymers, Polymer, 52, (2011),
5867-5878.
(13) Yashiro, K., Naito, M., Minagawa, Y. and Tomita, Y., On the Hysteresis of
Polyethylene and Polybutadiene: A Molecular Dynamics Study, Proceeding of
ICCES2005, Advances in Computational and Experimental Engineering and Science, CDROM(14) 屋代如月, 内藤正登, 皆川康久, 冨田佳宏, 非晶性高分子材料のヒステリシスに関する分子動力学的研究, 日本機械学会論文集, A72, (2006), 277-284.

(15) O’Brien, J. et al., An NMR Investigation of the Interaction between Carbon Blackand cis-Polybutadiene, Macromolecules, Vol.9, No.4(1976), pp.653-660.
(16) . Gert Heinrich, Manfred Klüppel, Thomas A. Vilgis., Reinforcement of elastomers, Current Opinion in Solid State and Materials Science, Vol.6, (2002),
pp.195-203.
(17) A. R. Payne., The Dynamic Properties of Carbon Black-Loaded Natural Rubber
Vulcanizates. Part I , Journal of Applied Polymer Science, Vol.6, (1962), pp.5763.
(18) J. A. C. Harwood, L. Mullins, and A. R. Payne., Stress Softening in Natural
Rubber Vulcanizates.Part II . Stress Softening Eﬀects in Pure Gum and Filler
Loaded Rubbers, Journal of Applied Polymer Science, Vol.9, (1965), pp.30113021.
(19) G. Kraus, C. W. Childers, K. W. Rollmann., Stress Softening in Carbon Black
Reinforced Vulcanizates. Strain Rate and Temperature Eﬀects, Journal of Applied
Polymer Science, Vol.39, (1966), pp.1530-1543.
(20) Dannenberg, E.M., The Eﬀects of Surface Chemical Interactions on the Properties
of Filler-Reinforced Rubbers, Rubber Chemistry and Technology, Vol.48, (1975),
pp.410-444.
63

(21) Julie Diani, Bruno Fayolle, Pierre Gilormini, A review on the Mullins eﬀect,
European Polymer Journal, Vol.45, (2009), pp.601-612.
(22) Julie Diani, Bruno Fayolle, Pierre Gilormini, On the Mullins eﬀect and hysteresis of ﬁbered biological materials: A comparison between continuous anddiscontinuous damage models, Estefania Pena, Juan A. Pena, Manuel Doblare,
International Journal of Solids and Structures, Vol.46, (2009), pp.1727-1735.
(23) F. Bueche., Molecular Basis for the Mullins Eﬀect, Journal of Applied Polymer
Science, Vol.4, (1960), pp.107-114.
(24) G. Marckmann, E. Verron, L. Gornet, G. Chagnon, P. Charrier, P. Fort., A theoryof network alteration for the Mullins eﬀect, Journal of the Mechanics and Physics
of Solids, Vol.50, (2002), pp.2011-2028.
(25) 内藤陽子, 伊藤眞義, シリカ充てん加硫ゴム材料の応力軟化に関する研究, 日本ゴム協会誌, 82-9,(2009),pp.394-399(26) 古谷泰大, 陸偉, 内藤正登, 冨田佳宏, 分子鎖網目モデルによるカーボンブラック
充填ゴムの粘弾性変形挙動シミュレーション, 日本機械学会第 17 回計算力学講演会講演論文集, No.04-40(2004-11), pp.17-18.

(27) 古谷泰大, 内藤正登, 陸偉, 冨田佳宏, カーボンブラック充填ゴムの繰り返し変形
挙動の評価, 日本機械学会論文集Ａ編, Vol.71, No.708(2005), pp.1109-1115.

(28) Tomita,Y.,Lu,W.and Furutani,Y., Micro- to Macroscopic Deformation Behavior of Carbon Black-Filled Rubber Under Monotonic and Cyclic Straining,
Proc.CIMTEC2004,PartB,(2004),pp.121-132.
(29) Tomita,Y.,Azuma,K.,and Naito,M.,Strain-Rate-Dependent Deformation Behavior of Carbon-Black-Filled Rubber under Monotonic and Cyclic Straining, KeyEngineering Materials,Vol.340,No.341(2007),pp.1017-1024.
(30) 児玉総治, 西勝志, 岩蕗仁, 古川睦久, シリカ配合天然ゴムの耐熱性, 日本ゴム協会
誌,Vol.76,No.11(2003),pp.393-398.

64

(31) Wolﬀ,S., Silanes in Tire Coupling After Ten Years - A Review, Tire Society andTechnology,Vol.15,No.4(1987),pp.276-294.
(32) 松沢憲治, シリカ/シランフィラーシステムの化学とゴム補強性, 日本ゴム協会
誌,Vol.78,No.6(2005),pp.211-217.

(33) 北村真瑠久, 神戸大学修士論文, シリカ充填ゴムの粘弾性挙動のモデル化と変形挙動の評価, (2010)(34) 望月利紀, 神戸大学修士論文, シリカ充填ゴムの内部微視構造のモデル化と変形応答のシミュレーションによる評価, (2011)(35) 中田伸哉, 神戸大学修士論文, FEM シミュレーションによるシリカ充填ゴムの微
視構造のモデル化と力学特性評価, (2013)(36) 田中文彦, 高分子の物理学, (1994), 裳華房.
(37) 冨田佳宏, ガラス状ポリマーの分子鎖網目理論による構成式と変形挙動のシミュ
レーション, 塑性と加工, Vol.37, No.424(1996), pp.485-491.

(38) Kuhn, W. and Grun, F., Beeziehuugen Zwischen Elastischen Konstanten und
Dehuungsdoppelbrechung Hochelastischer Stoﬀe, Kollooidzeitschrift, Vol.101,(1942), pp.248-271.
(39) 土井正男, 小貫明, 高分子物理・相転移ダイナミクス, (1992), 岩波書店.
(40) James, H.M. and Guth, E., Theory of the Elastic Properties of Rubber,J.Chem.Phys., 11 (1934), 455-481.
(41) James, H.M. and Guth, E., Theory of the elastic properties of rubber, Journal ofChemical Physics, Vol.11, (1943), pp.455-481.
(42) Wang, M.C. and Guth, E., Statistical Theory of Networks of Non-Gaussian Frexible Chains, Journal of Chemical Physics, Vol.20, No.7(1952), pp.1144-1157.
(43) Treloar, L.R.G., The Elasticity of a Network of Long-Chain Molecules-III, Transactions of the Faraday Society, Vol.42, (1946), pp.83-94.
65

(44) Flory, P.J. and Rehner, J., Statistical Mechanics of Cross-Linked Polymer Networks, Journal of Chemical Physics, Vol.11, No.11(1943), pp.512-520.
(45) Arruda, E.M. and Boyce, M.C., A Three-Dimensional Constitutive Model for
the Large Stretch Behavior of Rubber Elastic Materials, Journal of MechanicalPhysics and Solids, Vol.41, No.2(1993), pp.389-412.
(46) Green, A.E. and Zerna, W., Theoretical Elasticity, (1968), Oxford University
Press.
(47) Tomita,Y.,Azuma,K.and Naito,M.,Int.J.Mech.Sci., Vol50, No.5,(2008),pp856868.
(48) Truesdell, C. and Noll, W., The Non-Linear Field Theories of Mechanics,
Springer-Verlag, (1965).
(49) Boyce, M.C., Weber, G.G. and Parks, D.M., On the Kinematics of Finite Strain
Plasticity, International Journal of Plasticity, Vol.37, No.5(1989), pp.647-665.
(50) 北川浩, 弾・塑性力学−非線形解析のための基礎理論−, (1987), 裳華房.
(51) 久田俊明, 非線形有限要素法のためのテンソル解析の基礎, (1992), 丸善.
(52) Boyce, M.C., Parks, D.M. and Argon, A.S., Inelastic Deformation of Glassy Polymers. Part I :Rate Dependent Constitutive Model, Mechanics of Materials, Vol.7,(1988), pp.15-33.
(53) Doi, M. and Edwards, S.F., The Theory of Polymer Dynamics, (1986), Oxford
University Press.
(54) Gennes, P.G. de., Reptation of a Polymer Chain in the Presence of Fixed Obstacles, Journal of Chemical Physics, Vol.55, (1971), pp.572.
(55) Raha, S. and Bowden, P.B., Birefringence of Plastically Deformed Polymethyl
Methacrylate, Polymer, Vol.13, (1972)，pp.174-183.

66

(56) Botto,P.A., Duckket,R.A. and Ward,I.M., The Yield and Thermoelastic Properties of Oriented Polymethyl Methacrylate, Polymer, Vol.28, (1987)，pp.257-262.
(57) Tomita, Y. and Tanaka, S., Prediction of Deformation Behavior of Glassy Polymers Based on Molecular Chain Network Model,International Journal of Solids
and Structures, Vol.32, No.23(1995), pp.3423-3434.
(58) Tomita, Y., Adachi, T. and Tanaka, S., Modelling and Application of Constitutive Equation for Glassy Polymer Based on Nonaﬃne Network Theory, European
Journal of Mechanics, A/Solids, Vol.16, No.5(1997), pp.745-755.
(59) Higa, Y. and Tomita, Y., Computational Prediction of Mechanical Properties of
Nickel-based Superalloy with Gamma Prime Phase Precipitates, Proceedings of
ICM8,Advance Materials and Modeling of Mechanical Behavior, Vol.III, (1999),
pp.1061-1066.
(60) 比嘉吉一, 冨田佳宏, 粒子強化型複合材の均質化法による変形挙動のモデル化とシ
ミュレーション, 日本機械学会論文集Ａ編, Vol.66, No.648(2000), pp.1441-1446.

(61) 冨田佳宏, 数値弾塑性力学, (1990), 養賢堂.
(62) Guedes, J.M. and Kikuchi, N., Preprocessing and Postprocessing for Materials
Based on the Homogenization Method with Adaptive Finite Element Methods,
Computer Methods in Applied Mechanics and Engineering., Vol.83, No.2(1990),
pp.143-198.
(63) 内田真, 冨田佳宏, 非晶性高分子材料の微視組織の 不均一性が変形挙動に及ぼす影響, 日本機械学会論文集 A 編, 70-698, (2004), 1484-1490.

(64) 内藤正登, 東圭佑, 冨田佳宏, フィラー充填ゴムの粘弾性挙動のモデル化と数値シ
ミュレーション, 日本機械学会第 19 回計算力学講演会講演論文集, No.06-9(2006-

11), pp.397-398.

67

付録 A
非圧縮性ゴム粘弾性体の構成式の速度形式表示
式 (2.10) の速度形を導出する．

1
σi =
3

{

(
)}
(
)
√
√
λ
1
λ
λ2i
c
β
CαR Nα L−1 √+CβR Nβ L−1 √
λγ
λc
Nα
Nβ
{
(
)}
′
√
λcB
1
λ2i
R
−1
√
+
CαB NαB L
−p
3
λc
NαA

式 (A.1) を左 Cauchy-Green 変形テンソル Aij を用いて表すと，
(
)}
{ √
(
)
√
λ
1
λ
1
λ2i
c
β
CαR Nα L−1 √σij =+CβR Nβ L−1 √
3
λγ
λc
Nα
Nβ
{
(
)} 2 ′
√
1
λcB
λi
R
−1
√
+
CαB NαB L− pδij3λc
NαAとなる．続いて，下記の関係式∇σ= σ̇ − W σ + σW ,I˙1 = (trA)· = 2A · D,Ȧ = (D + W ) A + A (D − W )を用いることにより，式 (A.2) の速度形式表示は次式になる．
[{
(
)}
√
)
R
′
′
√ ( ζ
C
N
1
L
L
ζ
∇
β
β
σ ij =√CαR Nα √+−Aij Akl /Amm
−
3
λγ
Nα λc
λγ Nβ λc{]√ }√LCαR Nα L′ CβR Nβ++{δik Ajl + Aik δjl } ε̇kl
λc
λc

68

(A.1)

(A.2)

(
)
√
[{
( ′′
)}
√
CβR Nβ λ̇γ1ζλβ ζ ′
L′′
R
′
− 2√
L +√
Aij +
CαB NαB √−A′ij A′kl /A′mm
3
λγ 3AmmNαB λcB
Nβ
√
]
R
}
NαB {
L′′ CαB+δik A′jl + A′ik δjl (ε̇kl − ε̇pkl ) − ṗδij
(A.3)
λcB
∇

ただし，σ ij は Cauchy 応力の Jaumann 速度で，I1 は左 Cauchy-Green 変形テンソル

Aij の第 1 主不変量で，I1 = λ21 + λ22 + λ23 となる．D は変形速度テンソルで，W はス
ピンテンソルである．tr はトレースを表す．
∇

最後に，体積一定の条件下で，式 (A.3) に示す Cauchy 応力の Jaumann 速度 σ ij を∇Kirchhoﬀ 応力の Jaumann 速度 S ij で置き換えても本質的な差はないことと，変形速度テンソル D をひずみ速度テンソル ε̇ で置き換えることにより，非圧縮性ゴム弾性体の構成式の速度形式表示は次式になる．

[{
(
)}
√
)
R
′
′
√ ( ζ
C
N
1
L
ζ
L
β
β
√
−
CαR Nα √+−Aij Akl /Amm
S ij =3λγNα λc
λγ Nβ λc]{√ }√LCαR Nα L′ CβR Nβ+{δik Ajl + Aik δjl } ε̇kl
+
λc
λc
(
)
√
[{
)}
( ′′
′
√
CβR Nβ λ̇γ
λ
ζ
L′′
1
ζ
β
′
R
L +√
−
− 2√
Aij +
CαB NαB √A′ij A′kl /A′mm
3λλγ 3Amm
NαB
Nβ
cB
√
]
R
}
NαB {
L′′ CαB+δik A′jl + A′ik δjl (ε̇kl − ε̇pkl ) − ṗδij
(A.4)
λcB
∇

√
d −1
ここで， N は分子鎖の限界伸び比を表す．ζ = dx
L (x)

69

√
x=

2

I1
3N

= 1−β 2βcsch2 β である．

付録 B[ϕ],[B],[E],{ψ} の具体形

3(x3,y3 )Ae2

Ae1
P(x,y)Ae3

1(x1 ,y1 )2(x2,y2 )Fig.B.1 Triangle Element

要素は図 B.1 に示すような三角形 1 次要素を用いるので，形状マトリクス [ϕ] 及び{ψ} は次のようになる．
{T } = {ψ}T {θ} , {ψ}T = {ψ1 ψ2 ψ3 }T{ }vx= [ϕ]{δ̇}{v} =vy[ϕ] = [ [ϕ1 ][ϕ2 ][ϕ3 ] ]][ψi 0[ϕi ] =
0 ψi
70

(B.1)
(B.2)
(B.3)
(B.4)

ここで，ψi は形状関数である．図 B.1 のように節点 1,2,3 の座標をそれぞれ (x1 , y1 ),(x2 , y2 ),(x3 , y3 ) とし，要素内の任意の点 P の座標を (x, y) とすると，全体の面積 Ae ,Ae1 ,Ae2 ,Ae3 は次のように表せる．

1 x1 y1
2Ae = det  1 x2 y2 
1 x3 y3

1 x y
2Ae2 = det  1 x3 y3 
1 x1



1

2Ae1 = det 1
1
1
2Ae3 = det  1
1

,

,

y1

x
x2
x3
x
x1
x2


y
y2 
y3

y
y1 

(B.5)

(B.6)

y2

上式を用いて形状関数 ψ1 ,ψ2 ,ψ3 は次式のようになる．

ψ1 =

Ae1
Ae

,

ψ2 =

Ae2
Ae

,

ψ3 =

次にマトリクス [B],[E] は平面問題では次のようになる．
 ε̇xx {ε̇} =ε̇yy= [B]{δ̇}
2ε̇xy[B] = [ [B1 ][B2 ][B3 ] ]


ψi,x
0
[Bi ] =  0
ψi,y 

ψi,y ψi,x
vx,x 



v 
y,y
= [E]{δ̇}
{q} =

v 


 x,y vy,x[E] = [ [E1 ][E2 ][E3 ] ]


ψi,x
0
 0
ψi,y 
[Ei ] =  ψi,y0 
0 ψi,x

71

Ae3
Ae

(B.7)

(B.8)
(B.9)
(B.10)(B.11)(B.12)(B.13)

付録 C
関連発表論文・講演論文
学術講演 本馬剛徳 ，屋代如月，冨田佳宏，FEM シミュレーションによるシリカ充填ゴムのモデル化と力学特性評価，日本機械学会 第 26 回計算力学講演会 CMD2013, 佐賀大学, (2013.11)

 本馬剛徳，屋代如月，冨田佳宏，FEM シミュレーションによるシリカ充填ゴムのモデル化と力学特性評価，日本機械学会 第 27 回計算力学講演会 CMD2014, 岩
手大学, (2014.11)72

論文貼り付け用ページ

73

論文貼り付け用ページ

74

論文貼り付け用ページ75

論文貼り付け用ページ76

謝辞
本研究を遂行するにあたり，著者が大学在籍時より研究全般に渡り懇切丁寧な御教授，御指導を賜るとともに，貴重な研究発表の機会や素晴らしい研究環境を提供して頂いた屋代如月教授 (現岐阜大学) に深甚な感謝の意を表します．また，日頃から研究について多くの御助言，御教示を賜るとともに，本論文の作成にあたり丁寧かつ熱心なご指導を賜りました冨田佳宏教授（現同志社大学）に心より感謝致します．
本研究に対して貴重な御助言と御指導，有益な議論を賜りました田中克志教授に深く感謝致します．また，日々の研究活動を円滑に行えるよう便宜を図っていただきました片山雷太技術職員に感謝致します．
本研究において実験データを提供していただいた住友ゴム工業株式会社の諸氏に厚く御礼申し上げます．同じ研究テーマに取り組み，有益な討論と適切なご指摘を頂いた，中田伸哉氏（現住友電工ハードメタル株式会社) に心から感謝致します．
日々の研究活動や，大学院 1 年次の就職活動中に多くの有益なご意見を頂いた表面・界面工学研究室の諸先輩方に深く感謝いたします．また研究室配属以来，同輩として互いに切磋琢磨し，時には励ましあった伊藤貴彦氏，小出卓矢氏，津川悠太氏，菱谷卓也氏，藤麻成貴氏，毛利友宙氏，森太一氏に厚く御礼申し上げます．表面・界面工学研究室諸氏に今日に至るまで数々の御支援と御協力を頂きましたので，ここに記し感謝の意を表します．
最後に，6 年間の長きに渡る学生生活において，多大な支援と理解を頂いた両親に，敬意と感謝の意を表します．
平成 27 年 2 月
本馬77

剛徳