修士論文

FEM シミュレーションによるシリカ充填ゴムの微視構造のモデル化と力学特性評価

指導教員：屋代 如月

中田 伸哉

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

Master Thesis
Finite Element Modeling and Characterization of
Microstructure in Silica-Filled Rubber

February 2013
Division of Mechanical Engineering,
Graduate School of Engineering,
Kobe University, Kobe, JapanSHINYA NAKATA

要約
本研究では，シリカ充填ゴムの粒子分布形態や界面ゲル相の物性について，種々の実験結果に基づくモデル化を非アフィン分子鎖網目理論と均質化有限要素法の枠組みの中で試みるとともに，最大応力やヒステリシスロス増大のキーメカニズムについて基
礎的な検討を行った．
まずシリカ充填ゴムの分散構造を単純な直鎖状数珠繋ぎ構造でモデル化し，粒子の配置形態が応答に及ぼす影響を検討した．ゴム相には実験により得られた応力ースト
レッチ関係を再現する非アフィンモデルを，ゲル相には模擬実験をベースに作成した
アフィンモデルを用いている．その結果，粒子を繋ぎ止めるゲル相に変形が集中し，系の応力上昇の主要因となることが分かった．ゲル相アフィンモデルではヒステリシス
ロスが小さくなるため，非アフィンモデルをゲル相に導入し，かつ初期セグメント数

N s を小さく（硬いゲルを仮定）したところ，分子鎖の配向硬化が強められ，かつ除荷時の応力低下が顕著になりヒステリシスロスが増大することが示された．次に，実験で新たに報告された硫黄架橋時に生じる絡み点の不均一性を考慮した解析を行った．
その結果，局所的に硬い部分の存在が，カーボンブラック充填ゴムの検討で既に報告している，見かけの体積減少と同様の効果をもたらし最大応力が上昇すること，また軟らかい部分で顕著な非アフィン変形を生じることでヒステリシスロスが増大することが分かった．最後に，不均一性の効果を考慮したゴム相・ゲル相の応答を非アフィンモデルの構成式で再フィッティングし，45◦ 数珠繋ぎ構造の FEM シミュレーションを行った．その結果，ゲル相の分子鎖の配向硬化がさらに顕著になり不均一性をもたない最初のモデルの解析結果よりも最大応力・ヒステリシスロスが増大した．

Summary
Toward a new FEM model for dispersion morphology of silica particles and mechanical properties of gel interface in silica-filled rubber, we constructed various constitutive
models in the framework of nonaﬃne chain network theory and the homogenization
method and investigated the key mechanisms for maximam stress and hysteresis loss.
First, we examined the influence of dispersion morphology of silica particles on themacroscopic stress-strain responses, with a simple straight-chain bunching structure
where the interfacial gel phase connects silica particles. The result shows that theremarkable stress increase occurred due to the deformation concentration on the connecting gel phase. Aﬃne chain-network model for gel phase shows small hysteresis
loss than experimental observation, so that we introduced nonaﬃne model for gel
phase and discussed the influence of the initial segment number N s (smaller segment
indicate harder gel), then revealed that the strain hardening by chain orientation is
activated, resulting in drastic stress decrease in the unloading process and increase in
the hysteresis loss. Then, we discussed about the eﬀect of the heterogeneity of

chain

entanglement , which is recently reported by experiment, in the rubber matrix andgel phase. As a result, we found that the heterogeneity enhances the stress increasesince the deformable part decreases due to the presence of hard part, as previously
reported in the analysis of carbon black filled rubber. It also causes the remarkable risein the hysteresis loss since significant nonaﬃne deformation takes place at soft part.
Finally, we re-fitted the viscoelastic constitutive models for rubber and gel phases taking the eﬀect of the heterogeneity into account. FEM simulation for 45◦ bunching
structure demonstrated that maximum stress and hysteresis loss are further increased
due to the remarkable orientation hardening in the gel phase, compared to the resultof homogeneous gel model.

目次
第 1 章 緒論

1

第 2 章 基礎理論

4

2.1

2.2

ゴム単相の構成式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4

2.1.1

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

4

2.1.2

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

7

2.1.3

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

14

均質化法による微視組織のモデル化 . . . . . . . . . . . . . . . . . . . .

16

2.2.1

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

16

2.2.2

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

19

第 3 章 直鎖状数珠繋ぎ構造をもつシリカ充填ゴムの粘弾性変形挙動

23

3.1

ゴムマトリクス部のモデル化. . . . . . . . . . . . . . . . . . . . . . .

23

3.2

ゲル相のモデル化 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

24

3.2.1

ゲル相物性見積もり実験 . . . . . . . . . . . . . . . . . . . . . .

24

3.2.2

アフィンモデルによるゲル相のモデル化 . . . . . . . . . . . . .

25

シリカ粒子の配置形態の基礎的検討 . . . . . . . . . . . . . . . . . . . .

27

3.3.1

解析モデル . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

27

3.3.2

0◦ 数珠繋ぎ構造モデルの解析結果 . . . . . . . . . . . . . . . . .

29

3.3.3

30◦ 数珠繋ぎ構造モデルの解析結果 . . . . . . . . . . . . . . . .

33

3.3.4

45◦ 数珠繋ぎ構造モデルの解析結果 . . . . . . . . . . . . . . . .

35

ゲル相の非アフィンモデル導入 . . . . . . . . . . . . . . . . . . . . . .

37

3.4.1

非アフィンモデルによるゲル相のモデル化 . . . . . . . . . . . .

37

3.4.2

非アフィンゲル相モデルによる 45◦ 数珠繋ぎ構造の解析結果 . .

40

3.4.3

ゲル相の初期セグメント数 N s を変化させた場合 . . . . . . . . .

42

3.3

3.4

i

目次
第 4 章 分子鎖セグメントの不均一分布が変形挙動に及ぼす影響

4.1

4.2

4.3

5.2

45

ゴム単相の変形挙動評価 . . . . . . . . . . . . . . . . . . . . . . . . . .

45

4.1.1

解析モデル . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

45

4.1.2

ゴム単相モデルの解析結果 . . . . . . . . . . . . . . . . . . . . .

47

ゲル単相の変形挙動評価 . . . . . . . . . . . . . . . . . . . . . . . . . .

53

4.2.1

解析モデル . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

53

4.2.2

ゲル単相モデルの解析結果 . . . . . . . . . . . . . . . . . . . . .

54

ゲル単相の平均セグメント数 N0s を変化させた場合 . . . . . . . . . . .

59

4.3.1

解析モデル . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

59

4.3.2

N0s =6.4,4.4 としたゲル単相モデルの解析結果 . . . . . . . . . .

59

第 5 章 分子鎖の不均一性を反映させた構成式の構築

5.1

ii

63

セグメント不均一性を反映したゴム単相・ゲル単相の構成式化 . . . . .

63

5.1.1

ゴム単相の構成式化 . . . . . . . . . . . . . . . . . . . . . . . .

64

5.1.2

ゲル単相の構成式化 . . . . . . . . . . . . . . . . . . . . . . . .

66

不均一性を反映した構成式による 45◦ 数珠繋ぎ構造モデルの解析 . . . .

68

第 6 章 結論

71

参考文献

74

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

80

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

82

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

84

謝辞

97

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

(a)Tire.

(d)Molecular chain.

100nm(c)CB filled rubber.

Fig.1.1 Hierarchy of rubber(4) .

1

力学特性を示す．それらの力学特性の発現メカニズムの解明には高分子材料の微視的構造の観察が必要不可欠であり，柴山らは中性子小角散乱（SANS）装置を用いたゴム
内部の硫黄架橋の不均一性の観察 (5) を行っており，中島らは原子間力顕微鏡（AFM）
を用いた不均一構造の可視化に関する研究 (6)(7) を推進している．池田らは，変形中のゴム内部の構造変化を観察することに成功しており，分子鎖の伸長に伴う結晶化について言及している (8) ．さらに，分子動力学（MD）を用いた高分子材の分子レベルの変形挙動解析も盛んに行われており (9)(10)(11) ，屋代らは分子鎖の配向による応力上昇
ならびにヒステリシス発現のメカニズムについて新たな知見を報告している (12)(13) ．
充填ゴムのヒステリシスや粘弾性応答は，前述の補強効果により未充填ゴムに比べて顕著になることが確認されており (2)(14)(15) ，充填ゴムに引張変形を加えた時に生じる応力軟化（＝ヒステリシスロスの増大）の要因として，フィラー粒子の凝集体構造の変化 (16)(17)(18) や，ゴムの延伸に伴う高分子鎖の滑り (19) ，架橋構造の破壊 (20)(21) が古くから考えられてきた．最近では内藤らが，ゴム分子鎖のエネルギー変化を詳細に検討することにより，変形中の架橋構造・凝集体構造の破壊を裏付ける実験結果を報告している (22) ．このように充填ゴムはその内部で非常に複雑な力学特性を呈することが知られており，製品の設計や製造に充填ゴムの力学特性を最大限に生かすためには，
材料内部の微視的変形挙動の解明と，フィラー/ポリマー間の相互作用の適切な評価，
そしてそれらを反映した巨視的応答を精密に再現し得るシミュレーションモデルの構築が求められている．我々の研究グループではこれまで，CB 充填ゴムの高機能性発現の詳細なメカニズム，特に CB 充填に伴うゴム部の微視的変形挙動と CB 充填ゴムの巨視的応答の関係を明らかにするため，ゴム部の適切な構成式の定式化と，CB 充填ゴムが材料の機械的特性に及ぼす影響，微視領域における変形挙動を検討し得るシ
ミュレーションモデルの構築，評価など，広範多岐に及ぶ研究を推進し，得られた成果はタイヤの実際の設計にも用いられている (23)(24)(25)(26) ．
現在，工業的に利用されているゴム材料の多くは多量生産が容易な CB が充填剤として用いられているが，その一方で石油を原材料とせず脱石油に貢献するシリカが新たな充填剤として脚光を浴びている．シリカ充填ゴムは CB 充填ゴムよりも耐油性・耐酸性に優れ，タイヤに用いる場合では転がり抵抗が小さくなるため燃費の向上をもた
らし高性能なタイヤ材料となり得ることが報告されている (27) ．最近の実験から，シリカ充填ゴムの製造工程においてシリカ粒子とゴム材の界面状態を制御するために添加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-filled rubber.
するシランカップリング剤の影響 (28)（図 1.2(b)）により，界面にゴムとは物性の異なる「ゲル相」が生成し，図 1.2(a) のようにシリカ粒子を数珠繋ぎ状に繋ぎ止めていることが観察された (4) ．そこで，我々の研究グループでは CB を対象として開発してきたこれまでのフィラー充填ゴムの FEM モデルを発展させ，シリカ粒子のネットワーク
構造や数珠繋ぎ構造，ゲル相の厚さなどを様々に変化させたシミュレーション (29)(30)を行ってきたが，実験結果に見られる引張後期に見られる配向硬化がもたらす応力上昇ならびにヒステリシスロスが十分に再現できていなかった．本研究では，数珠繋ぎ構造を呈するシリカ充填ゴムの力学応答の基礎を把握するものとして，シリカ粒子が直列に連結した単純な解析モデルを用いて，最近の実験により報告された硫黄架橋の不均一性 (5) を考慮したモデル化を行った．
第 2 章では解析手法の基礎として，ゴムの粘弾性応答を表現するための分子鎖網目
理論に基づく粘弾性 8 鎖モデルの構成式について説明し，解析に用いる有限要素均質化方程式の概要を述べる．第 3 章では，実験データを基に同定したゴムマトリクス部
およびゲル相の粘弾性 8 鎖モデルの構成式を導入した直鎖状数珠繋ぎ構造の FEM シ
ミュレーションを行い，シリカ粒子の配置形態，ならびにゲル相の特性変化が系の応答に及ぼす影響を明らかにする．第 4 章では，硫黄架橋により生じる絡み点の空間的な不均一性を，分子鎖網目理論におけるセグメント数 N の不均一性として導入し，分布形態の違いがゴム単相，ゲル単相の応答に及ぼす影響を検討する．第 5 章では，第4 章の解析結果を粘弾性 8 鎖モデルの構成式として再度パラメータフィッティングする
ことでスケールアップを行い，第 3 章の解析結果と比較し議論する．最後に，第 6 章で本研究で得られた結果の総括を述べる．

3

第2章
基礎理論
本章では，まずゴム弾性応答を記述するために提案された分子鎖網目理論 (31)∼(40) について説明する．次に，分子鎖に管模型 (48) を用い，非ガウス鎖理論に基づく非圧縮性を考慮したゴム粘弾性体の構成式について説明する．さらに，シリカ充填ゴムにおける分子鎖のからみ点数の変化を許容する非アフィンモデルへの一般化を行う．次に，

2 変数漸近展開理論に基づく均質化法の基本的な考え方を述べた後，ゴムの構成式を更新ラグランジュ法に基づく均質化理論 (54, 55) に導入することにより，微視的関係式
及び巨視的平衡式を導出し，その有限要素表示式を示す．

2.1

ゴム単相の構成式2.1.1

分子鎖網目理論

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

θ を一定に保ちながら，第 2 結合を軸に回転可能となり，その位置は，回転角の関数として表されるポテンシャルエネルギによって決まる．このような考え方で高分子材料

4Fig.2.1 Concept of hierarchical structure(a)molecular chain, (b)segment,(d)macroscopic continuum.

of polymeric material(32) ,(c)structure of network,の微視的構造を忠実に考慮したモデルを構築し，高分子材料の挙動を表現することが原理的には可能である．しかしながら，実際の適用に当たっては多くの時空的な制約が加わるため，モノマーを直接扱わず，セグメントという最小構成単位で粗視化した分子鎖網目モデル (32) が用いられる．
分子鎖網目理論では，高分子材料は図 2.1(c) に模式的に示すように，分子間の化学的結合あるいは物理的結合により接合点において連結された鎖が，ランダムに配向した網目構造を有していると仮定している．さらに，(i) 接合点は原子の揺らぎ周期に対して長時間的には平均位置が変化せず，接合点周りの摂動は無視できる，(ii) 二つの接合点を両端に持つ分子鎖の端−端ベクトル (end-to-end vector) は，それが埋め込まれている材料の連続体と共変形をするとの仮定を置く．このようなモデルをアフィンモ
デルという．図 2.1(b) に示すように，二つの接合点間の分子鎖は「1 本の分子鎖」と定義され，それは導入した最小構成単位である「セグメント」から成る．モノマーの数が十分多ければスケーリング則によって鎖の巨視的な性質は変わらない (34) ．ここで導入した「セグメント」が「モノマー」とどう対応するかは議論の対象とせず，ここでは現象論的に最小単位のセグメントを導入している．
「1 本の分子鎖」の形態が非ガ
ウス統計分布 (35) に従うとすると，二つの接合点を結ぶ方向にストレッチ λ を加えた

5

Fig.2.2 Molecular chain network model, (a) three chain model, (b) four chain
model, (c) eight chain model.

場合に生じる応力 σ は次式で表すことができる (33) ．
(
)
√
λ
−1
√
σ = kB T N L
N

(2.1)

ここで，N は 1 分子鎖あたりのセグメント数，kB は Boltzmann 定数，T は絶対温度である．また，関数 L(x) は次式で定義される Langevin 関数である．
{ ()}
sinh xd
1
L(x) =ln= coth x −
dx
x
x

(2.2)

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

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

6

2.1.2

ゴム粘弾性体の構成式

図 2.2(c) で示される８鎖モデルはゴム超弾性体の変形応答を記述するのに用いられている．超弾性体とは負荷を受け大きく変形した後，完全に除荷すると元の状態に戻る弾性体である (41) ．しかしながら，実際の分子鎖は周囲の分子鎖からの摩擦に起因する粘性も持ち合わせている．そこで，ゴムの粘弾性挙動を記述するために，図 2.3 に示すような粘弾性８鎖モデルとダンパーで構成されるモデルを構築する．

B
A

β
α
γ

D

Fig.2.3 Revised eight chain model.

周囲の分子鎖からの摩擦を表現するために，図 2.2(c) に示す８鎖モデルの各単鎖に，
図 2.3 挿入図に示す粘性抵抗をもつバネ・ダンパーの標準モデルを導入した新たな 8 鎖
モデル A が提案されている (42) ．ここで用いた粘性抵抗は，後に説明する管模型 (48) によって表される．ここで，初めに８鎖モデル A の構成式を記述する．式 (2.1) より，単分子鎖の二つの接合点を結ぶ方向にストレッチ λc を加えた場合，図 2.3 挿入図に示す
システムにストレッチ λc を加えた場合に生じる応力 σc は次のように表せる．

σc = CαR√(Nα λc L−1

λ
√c
Nα

)

√
R

+ Cβ

λc
Nβ L−1
λγ

(

λ
√β
Nβ

)

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

(2.3)
√
√
Nα ， Nβは分子鎖の限界伸びを表す．nα ，nβ は構成要素のバネ α，β に含まれる分子鎖の数を表す．また，図 2.3 に示す単分子鎖の各要素のストレッチを λα ，λβ ，λγ とし，λα = λcである．その他の添え字 α，β ，γ についても，図 2.3 に示す要素 α，β ，γ と対応して7

いる．ただし，λc = λβ λγ(2.4)

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

σ c = λc

∂W
∂λc

(2.5)

式 (2.3)，式 (2.5) より恒等的に次式が成り立つ．

√∂Wc= CαR Nα L−1
∂λc

(

λ
√c
Nα

)

√
1
+ CβR Nβ L−1
λγ

(

λ
√β
Nβ

)
(2.6)

8 鎖モデル (40) の場合，主ストレッチを λ1 ，λ2 ，λ3 とすると，分子鎖のストレッチは√
λc = (λ21 + λ22 + λ23 ) /3 と表すことができるので，
∂λc
λi
=
∂λi
3λc

(2.7)

の関係が成り立つ．式 (2.5)-(2.7) より，8 鎖モデル A の主ストレッチ方向の応力 σiA とストレッチ λi は次の関係で与えられる．

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

σiA = λi(2.8)

ここで，図 2.3 の８鎖モデル B の主ストレッチ方向の応力 σiB と弾性ストレッチ λ′i は次のように表せる (40) .

1
σiB =

3

{
R
CαB

√

(
−1

NαB L

λ
√ cB
NαB

)}

λ2i
λcB

′

(2.9)

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

ゴム粘弾性体の変形は体積変化が小さいとしてそれを無視する場合が多い．そこで本
研究では，非圧縮性ゴム粘弾性体を取り扱うものとし，非圧縮性を満たすために静水
圧 p を用いる．この時，式 (2.8) , (2.9) を用いると，非圧縮性ゴム粘弾性体の構成式は次式のように表せる．

σi = σiA + σiB − p

8(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β λ̇γ
λ
ζ
1
ζ
L′′
β
′
R
− 2√
L +√
Aij +
CαB NαB √−A′ij A′kl /A′mm
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を次式で定義する (44) ．

F =

∂x
∂X

(2.12)

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

F = F βF γ

(2.13)F γ は完全な除荷状態で応力解放配置を表す．また，変形勾配 F は弾性ストレッチ V β ，
回転 R，粘性ストレッチ U γ を用いて次の形で表現される．

F = V β RU γ

9(2.14)

F
Current
Configuration
Reference
Configuration

F

β-1

Fβ

RelaxedConfiguration
γ

F

Fig.2.4 Concept of viscoelasitc decomposition of deformation gradient.

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

R = Rβ Rγ ,(2.15)

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

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)とすることによって，次式の関係を得る．

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

L=∂v= d + w = Ḟ F −1 = Ḟ β F β−1 + F β Ḟ γ F γ−1 F β−1∂x10(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)σ γ は要素 γ にかかる応力である．式 (3.1) より，弾性ストレッチ λβ は分子鎖一本のス
トレッチ λc を用いて，次のように表される．

λβ =

λc
λγ

粘性ストレッチ λγ は次のように表現される．
∫ tλγ =λ̇γ dt(2.26)(2.27)0˜ γ|
λ̇γ |t=t+∆t = λγ |t=t ·|d(2.28)

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

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

σ ∗γ = σ γ
11(2.30)

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

σ ∗′
γ
N=√
2 τ∗

(2.31)

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

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

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

Bergström·Boyce は reptation 理論を基に，与えられた有効せん断応力 τ ∗ に対する

aaFig.2.5 Tube model.

12Undeformed Network

Deformed Network

Deformed and Relaxed Network

Fig.2.6 Relaxation behavior of polymer chain.

粘性せん断ひずみ速度 γ̇ γ を次のように導出した (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
ここで，Ĉ1 ，C2 ，m は材料定数であり，一般にひずみ速度に依存する．

13(2.37)

Physical Linkage

Initial State of Molecular Chain

Deformed State of Molecular Chain

Fig.2.7 Concept of deformation of molecular chain accompanies decrease inphysical linkage.

2.1.3

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

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

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

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

nα =

CαR
,
kB T

nβ =

CβR
,
kB T

14

nB =

R
CαB
kB T

(2.40)

Junction point

1

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

とする．ここで，Nn は系全体の総セグメント数，N0 は初期セグメント数であり，総セグメント数 Nn は材料に固有の数値である．また f (λc ) は λc の 2 次多項式で表し，f (λc ) = a0 + a1 λc + a2 λ2c(2.41)とする．a0 ，a1 ，a2 は定数である．
次にその f (λc ) の具体的な関数の決定方法について説明する．まず第一回目のサイクルでは，負荷時のからみ点数が変化し，除荷時はそれが変化しないとして，式 (2.10)を用いたシミュレーション結果と実験結果の差が小さくなるように，λc 関数の係数を同定する．一旦除荷した後，再負荷時においては，からみ点数の変化は不可逆なものとし，本シミュレーションでは再負荷時において前回のサイクルで到達した最大スト
レッチより小さい変形領域では平均セグメント数 N は変化しないものとした．すなわ
ち 1 回目のサイクルの除荷と 2 回目のサイクルの再負荷の N の値は同一とした．さらに変形が進み，前回のサイクルで経験した最大ストレッチを超えると，再び式 (2.38)
により N が変化するものとする．

15

2.2

均質化法による微視組織のモデル化

2.2.1

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

本節では，2 変数漸近展開理論に基づく均質化法の基本的な考え方を簡単に述べた
のち，ゴム材の構成式を更新ラグランジュ法に基づく均質化法 (54, 55) に導入することにより，微視的関係式及び巨視的平衡式を導出する．
図 2.9 に示すような全体構造 X の任意点の近傍において，局所的に周期性をもつ
微視構造 Y が存在する材料を仮定し，構造物全体を表現する座標系 xi (i = 1, 2, 3) とyi = xi /η の関係を満足する微視構造を表現する座標系 yi の 2 変数を導入する．ここで，η は微視的周期構造内の基本単位領域のスケールを表す．現変形状態における物体の体積を Ω，表面積を S ，外部表面の一部 St 上に作用する表面力を P とし，残りの外部表面 Su に一定の変位速度を与える．このとき，更新ラグランジュ法を用いると，
仮想仕事原理式は下記のように表せる (56) ．

∫ (

∫

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

St

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

Ṡij = S ij −Fijkl ε̇kl ,

W
xi

P

X

unit cellStYSu

yi = xi /hyi

Fig.2.9 The relation between overall structure and microscopic structure.

16

Fijkl =1(σlj δki + σkj δli + σli δkj + σki δlj )2(2.43)を用いることによって，次のような形に統一的に示すことができる．

Ṡij = Lijkl ε̇kl

(2.44)

一方，微視領域内の任意点の変位速度 v はスケールパラメータ η により，次のように漸近展開できる．

v = v (x, y) = v 0 (x, y) + ηv 1 (x, y) + η 2 v 2 (x, y) · · ·(2.45)式 (2.45) をひずみ速度と変位速度の関係式

1ε̇ij =
2

(

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

[]
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

(2.47)

次に，式 (2.44)，(2.45)，(2.47) を式 (2.42) に代入し，η について同じ次数の項を整理すると，以下の式が得られる．

1
η2

1
η

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

∫ {

[(

) ∂δv

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

∫ {
Ω

[(
Lijkl
[(

+σmj

) ∂δv

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

i

(
+

]

i

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

dV = 0) ∂δv ]

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

i

∂xj
∂yj
)
( 1
)
]}
∫
0
1
2
∂vi
∂vi ∂δvi∂vi∂vi ∂δvi+++dV =Ṗi δvi dS∂xm ∂ym ∂xj∂xm ∂ym ∂yj
St
17(2.49)(2.50)

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

( )
∫ ∫
x
1
dΩ →
Ψ
Ψ (y) dY dΩ,η|Y | Ω YΩ(2.51)( )∫ ∫1x
Ψ (y) dSdΩlim η ΨdS →η→0+
η|Y | Ω S
S(2.52)

∫
lim

η→0+

∫

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

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

(2.53)

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

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

∂v 0
∂v 0
Lijkl k + σmj i∂yl∂ym(2.54))
nj = 0(2.55)

Guedes ら (57) の命題 1 に基づき，式 (2.54)，(2.55) から，次式を得る．

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

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

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

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

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

]

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

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

2.2.2

有限要素均質化方程式

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

∫

(

T

T

∫

)

(

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

ここで L = D − F である．各マトリックスは次のように表せる．
)
(
)T
(
31
23
12
33
22
,
ε = ε11 ε22 ε33 2ε12 2ε23 2ε31
, χ,y = χ11
χ
χ
χ
χ
χ
,y
,y
,y
,y
,y
,y

χij
,y =

(

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

19

)T
,

1
χkl
(ij) =
2

(

∂χkl
∂χkl
j
i
+
∂yj
∂yi

)
,

(
q=

)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) =χij,y(q) =(

σ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
σyy

0
1
0
0
0
0

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)

χ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

sym.

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

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

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

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

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

T
ḋ1

T
ḋ2

···

,

T
ḋN =

T

χij
(d)N =()()

χ11 χ22 χ33 χ12 χ23 χ31

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

T

)T

T
ḋN

, χij
(d) =

(

, χij =
χij
(d)1

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

d˙N1 d˙N2 d˙N3

T

(

,

χij
χij
χij
1
2
3

χij
(d)2

T

)
,

· · · χij
(d)N

T

)T

)
.

T

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

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

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

δ ḋ

T

[∫

(

∫

)

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

T

Y

(

T

T

)

B L + E QR dY]=0(2.63)Y

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

[∫

(

T

T

]

)

∫

B LB + E QE dY χ(d) =( T)
B L + E T QR dY(2.64)YY

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

∫

(

T

H

T

H

)

δε L ε + δq τ q dV =
Ω
∫
(
)
1
H
L − Lχ,y dV,
L =|Y | Y
∫
(
)
1
H
τ =
Q − Qχ,y(g) dV|Y | Y21∫
δv T Ṗ dS,
St(2.65),

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

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

,

32
31
23
21
13
12
33
22
χ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 =

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

T

[∫

(

)

H

T

H

]

∫

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

Y

T

Ψ Ṗ dS = 0,St

(2.66)

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

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

(2.67)

St

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

22

第3章
直鎖状数珠繋ぎ構造をもつシリカ充
填ゴムの粘弾性変形挙動
3.1

ゴムマトリクス部のモデル化

True Stress σ22[MPa]6Experiment
Simulation

4

2

0

1

2

3

4

Stretch λFig.3.1 Comparison of true stress-stretch relations by simulation and experiment.

未充填ゴムにおける実験結果，ならびにそれを再現するようにパラメータフィッティ
ングしたシミュレーションでの結果を図 3.1 に示す．図 3.1 において，実験結果は変形
速度 u̇ = 100[mm/min] で最大ストレッチが λ = 4.0 になるまで 2 回繰り返し変形を与えた未充填ゴムのデータである．実際のゴムの変形応答にみられる除荷時の応力軟化
挙動や再負荷時の応力回復挙動を再現するため，シミュレーションは実験と同条件の

23

2 回繰り返し変形で，2.3 節で説明したからみ点数変化を許容する非アフィン変形を 8鎖モデル A 内の要素 α のみに適用し，要素 β や，８鎖モデル B のからみ点数は変化しないものとした．フィッティングにより得られた材料定数を表 3.1 に示す．構成式中の総セグメント数 Nｎ (= nα Nα + nβ Nβ + nB NB ) は Nｎ = 1.95 × 1027 である．また，材料の温度は変形過程を通して一定で 296[K] とした．

Table 3.1 Parameter of rivised eight chain model (Rubber matrix).

Rubber matrix
CαR (= nα kB T )CβR = (nβ kB T )

0.219
Ĉ1A
5.0 × 105

0.25

CBR (= nB kB T )
0.10

Nα
14.0

Nβ
14.0

NB
14.0

C2A
−0.50

mA
3.2

Ĉ1D
3.0 × 105

C2D
−0.50

mA
4.8

3.2

ゲル相のモデル化

3.2.1

ゲル相物性見積もり実験
mesh cage

silica filled rubber
gel and silica

toluenepress
rubber

Fig.3.2 Overview of experimental procedure of evaluation of gel phase.

これまでゲル相の物性は実験的に明らかにされていないが，カップリング剤の影響によりからみ点数が増加していると示唆されていることから，本研究ではゲル相における 1 分子鎖当たりのセグメント数 N が未充填ゴムに比べ少なくなっていると仮定し，それ以外の材料定数については未充填ゴムと同等として解析を行ってきた．ここでは，

24

True Stress σ[MPa]Gel phase
Unfilled rubber

2

1

0

1

1.2

1.4

1.6

1.8

2

Stretch λFig.3.3 Comparison of true stress-stretch relations for gel phase and unfilled rubber.

住友ゴム工業株式会社の新たな実験により見積もられたゲル相の物性を導入する．本項ではその実験手順について簡単に説明する．
図 3.2 に実験手順の概要を示す．まず初めにカップリング剤とシリカ粒子のみを未充填ゴムに混ぜ合わせ作成したシリカ充填ゴムを，ゴムの良溶媒であるトルエンに浸せきさせ，不溶分として残ったゲル相とシリカ粒子のみを抽出する．そして，取り出した粒子とゲルをプレスした後，温度分散を測定することによりゲルのガラス転移温
度 Tg を測定する．ガラス転移温度 Tg はゴム中の硫黄量に応じて変化する特性がある
ことから，測定したガラス転移温度 Tg の移動量により，ゲル相がどの程度の架橋密度に相当するかについて見積もりを行った．そして，見積もった硫黄量を添加し作成したゴムについて引張実験を行った．図 3.3 に見積もりを行った硫黄量を添加し作成したゴムと未充填ゴムに変形速度 u̇ = 100[mm/min] で最大ストレッチが λ = 2.0 になるまで 2 回繰り返し変形を与えた時の真応力ーストレッチ関係を示す．これより，ゲル相は未充填ゴムに比べ硫黄量が増加しているため硬化していることがわかる．

3.2.2

アフィンモデルによるゲル相のモデル化

3.2.1 項のゲル相物性見積もり実験では λ = 2.0 以上の試験は難しくデータが取れな
いが，シリカ充填ゴムの内部ではシリカ粒子を繋ぎとめるゲル相に大きな変形が集中していることが示唆されており，λ = 2.0 以降の挙動も非常に重要であると考えられ

25

True Stress σ22[MPa]25

Experiment
Gel phase 1
Gel phase 2
Gel phase 3

20

15

10

5

0

1

2

3

4

Stretch λFig.3.4 Comparison of true stress-stretch relations of three gel phase models
and experiment.

る．そこで本節では変形初期においてはこの実験結果にできるだけ従うようにしつつ，変形後期のゲル相の応力上昇挙動が異なる 3 種類のモデル化を行った．各モデルのパラメータを表 3.2 に示す．図 3.4 はゲル相物性見積もり実験結果と，作成した 3 つのゲル相モデルでの λ = 4.0 までの応答である．図のように実験により見積もられるヒス
テリシスロスは極めて小さいため，変形過程には 8 鎖モデル A，B ともに変形中の分子鎖のセグメント数変化を許容しないアフィン変形を適用している．
先述のようにゲル相はゴムマトリクス部に比べ硬い相となっていることが示唆されており，ゴム（ゲル）が硬い物性であるということは分子鎖同士がより複雑に絡まりあい分子鎖の絡み点数が多い状態であると考えられる．2.1.1 節で説明した分子鎖網目
理論に基づくと，分子鎖の絡み点数が多いということはすなわち 1 分子鎖当たりの平均セグメント数 N が小さいことに対応する．よってゲル相のセグメント数 N s はゴム
マトリクス部のセグメント数 N = 14.0 に対して N s = 8.4 とした．またゲル相はシリカ充填ゴム作成時に添加するカップリング剤の影響によってゴムが特性変化した物質であると考えられるが，ゲル相を構成する内部の分子鎖自体に変化はなく系全体の分子鎖の総セグメント数はゴム部と同じであるということを考慮し，ゲル相の総セグメ
ント数 Nｎ (= nα Nα + nβ Nβ + nB NB ) はゴムマトリクス部と同じ Nｎ = 1.95 × 1027 としてモデル化を行った．

26

Table 3.2 Parameters of rivised eight chain model (Gel phase, Aﬃne model).

Gel phase 1
CαR (= nα kB T )0.6483CβR = (nβ kB T )0.20

CBR (= nB kB T )
0.10

Nα
8.4

Nβ
8.4

NB
8.4

Ĉ1A
5.0 × 105

C2A
−0.50

mA
3.2

Ĉ1D
3.0 × 105

C2D
−0.50

mA
4.8

CαR (= nα kB T )CβR = (nβ kB T )CBR (= nB kB T )0.4283

0.42

0.10

Nα
8.4

Nβ
8.4

NB
8.4

Ĉ1A
5.0 × 105

C2A
−0.50

mA
3.2

Ĉ1D
3.0 × 105

C2D
−0.50

mA
4.8

Gel phase 2

Gel phase 3
CαR (= nα kB T )0.2183CβR = (nβ kB T )0.65

CBR (= nB kB T )
0.10

Nα
8.4

Nβ
8.4

NB
8.4

Ĉ1A
5.0 × 105

C2A
−0.50

mA
3.2

Ĉ1D
3.0 × 105

C2D
−0.50

mA
4.8

3.3

シリカ粒子の配置形態の基礎的検討

3.3.1

解析モデル

実際のシリカ充填ゴムは，シリカ粒子をランダムに分散した状態で含有するが，本節ではシリカ粒子の配置がシリカ充填ゴムの機械的特性に及ぼす基本的な影響を検討するために，平面ひずみ条件下でシリカ粒子が直鎖状に繋がった周期構造を対象とす
る．均質化法におけるミクロ構造の単位ユニットセル中に，2 個のシリカ粒子が界面ゲル相によって連絡された直鎖状の構造を有する解析モデルを図 3.5 に示す．図では界面ゲル相を水色で示している．直鎖状に配置したシリカ粒子の周囲の界面ゲル相の力学特性，ならびに粒子の配置が全体の変形挙動に与える影響のみを抽出するため，全てのモデルにおいてユニットセル内のシリカ粒子径は全て同じとし，シリカ粒子含有率は 30 ％，ゲル相割合は 12.54 ％ に統一し，2 つのシリカ粒子の配置を 0◦ から 45◦ まで変化させた 3 つの解析モデルを作成した．そのため，直鎖状 (0◦ ) からの傾きが大き

27

0c
L0

x2

y2

x1

0

y1

0

Unit cellx2L0(a) 0cbunching structure model

30c

y2

x1

0

0

y1

Unit cellx2L0(b) 30cbunching structure model

45c
y2

0

x1

0

y1

Unit cell(c) 45cbunching structure model
Fig.3.5 Simulation models of silica filled rubber.

くなるにつれ，ユニットセルも横に広がった形状となっている．
解析には，周期的微視構造を有する材料全体を表現する座標系 xi と微視構造を表現する座標系 yi の二変数を用い，変位を漸近展開する均質化理論 (57) に基づき定式化し28

た有限要素法 (54) を用いる．本研究で直接用いる均質化理論，並びに，有限要素方程式の具体形，計算手順の詳細については 2.4，2.5，文献 [5,15] を参照されたい．
巨視的に一様な単軸変形を与えるものとして，巨視座標系の x2 方向に u̇ = 100[mm/min]の一定な変形速度で，最大伸びが λ2 = 1.5 になるまで引張変形を与えた後除荷を行い，1 サイクルの変形挙動を解析する．ゴムの非圧縮性を満足させるためのペナルティ定
数は ṗ = 100 とする．シリカ粒子の剛性は，ゴム材の剛性に比べて十分大きいと考え，
ゴムマトリクス部とゲル相の計算の安定性と結果にほとんど影響を与えない値として，
縦弾性係数 E = 100[MPa], ポアソン比 ν = 0.3 とした．材料の温度は変形過程を通して一定で T = 296[K] とした．

3.3.2

0◦ 数珠繋ぎ構造モデルの解析結果

Nominal Stress Σn22[MPa]3Experimentwith Gel phase 1
with Gel phase 2
with Gel phase 3

2

1

0

1

1.2

1.4

Stretch λ2
Fig.3.6 Comparison of nominal stress-stretch relations of silica-filled rubber
by simulation and experiment (0◦ bunching structure model).

図 3.6 に，0◦ モデルにおいて界面ゲル相の物性を Gel phase 1 から Gel phase 3 まで 3パターンに変えて解析を行った結果を，実験から得られたシリカ充填ゴムの公称応力ストレッチ関係と比較して示す．今回実験におけるカップリング剤含有率 µ は，最も一般的に用いられる µ = 8[wt ％] とした．ここで，実験データはシリカ含有率 20 ％ の29

シリカ充填ゴムの引張試験結果であるのに対し，本章で用いるシリカ充填ゴムの解析
モデルのシリカ含有率は 30 ％ と，シミュレーションの方が高い応力を示しやすい条件
であることに注意されたい．いずれのモデルにおいても応力が大きく上昇し，λ2 = 1.5に達する前に計算不能となった．本モデルのようにシリカ粒子が縦に連結した条件では，界面ゲル相に非常に大きな変形が集中するため変形後期ではゲル相のメッシュが潰れ，計算不能となる．計算停止直前の λ2 = 1.43 までの解析結果を実験データと比較すると，当然ではあるが粒子を直列に連結して界面ゲル相の硬化が顕著に表れるようにした本モデルではいずれも実験値より高い値を示している．また本モデルでは，これまでのシリカ充填ゴムの界面相モデル (29) ，ネットワークモデル (30) では再現できなかった，実験結果にみられる変形後期（λ2 = 1.35〜1.50）での応力上昇挙動（3 次関
数的な負荷曲線）の傾向を再現可能となっている．

λc
3.30
2.84
2.38
1.92
1.46
1.00

(a)Gel phase1(b)Gel phase2(c)Gel phase3

Fig.3.7 Distribution of morecular chain stretch λc (0◦ bunching structure model).
計算停止直前の λ2 = 1.43 における分子鎖ストレッチ λc の分布を図 3.7 に，引張方
向応力 σ22 の分布を図 3.8 に示す．0◦ モデルにおいてはシリカ粒子が直線的に連結されており，剛性の高いシリカ粒子はほとんど変形しないためユニットセルモデルの中心軸付近では系全体に与えた引張変形がゲル相に集中し，シリカ粒子同士を縦に連結しているゲル相において分子鎖ストレッチ λc が非常に大きくなっている．またユニット
セルモデルの周期性のため 2 つの粒子に挟まれた粒子の側面部分のゴムマトリクス部
では系全体の延伸に伴い粒子から圧縮変形を受けるが，ゴムは体積一定を満たしながら変形するためその分縦に伸び，他のゴム部分に比べると大きな分子鎖ストレッチ λc30

σ22[MPa]
58.40
46.69
34.98
23.27
11.56
-0.15

(a)Gel phase1(b)Gel phase2(c)Gel phase3

Fig.3.8 Distribution of tensile stress σ22 (0◦ bunching structure model).

を生じている．また 8 鎖モデルの幾何学関係から，図 3.4 の横軸に示したゲル相モデ
ル全体の引張変形 λ と λc の間には

√ (
)
1 1
2
+1+λ
λc =
3 λ2

(3.1)

の関係が成り立つ．したがって図 3.7 において最も分子鎖ストレッチの大きい部分では

λc = 3.30 となっているが，これをゲル相モデル全体の引張変形 λ に換算すると λ = 5.68に相当する．このためシリカ粒子を縦に連結しているゲル相には図 3.4 に示した以上の引張変形が加わっており，図 3.8 のように非常に大きな応力を生じる．

λ2 = 1.43 までの引張方向応力 σ22 分布の変化を図 3.9 に示す．いずれのゲル相物性を導入したモデルにおいても，(a) → (c) の間は粒子同士を繋ぎ止めているゲル相の応力上昇が比較的なだらかであるのに対し，(c) → (e) にかけては急激に上昇している．
これは，ゲル相に分子鎖ストレッチ λc の増大とともにセグメント数 N s の変化を生じないアフィン変形を適用しているためであり，変形後期に非常に強い配向硬化を示す．
このようなシリカ粒子間のゲル相で生じる顕著な応力上昇が，図 3.6 で示した変形後期での 3 次関数的な応力上昇挙動をもたらしたと考えられる．またいずれの場合も粒子を数珠繋ぎ状に繋ぎ止めていたゲル相が系の延伸とともに細長く引き伸ばされ，ある程度の距離をあけて粒子を繋ぎ止めるネットワーク形状となっている．

31

σ22[MPa]
58.40
46.69
34.98
23.27
11.56
-0.15
(a)λ=1.10(b)λ=1.20(c)λ=1.30(d)λ=1.40(e)λ=1.43(d)λ=1.40(e)λ=1.43(d)λ=1.40(e)λ=1.43

Gel phase 1

σ22[MPa]
58.40
46.69
34.98
23.27
11.56
-0.15

(a)λ=1.10(b)λ=1.20(c)λ=1.30

Gel phase 2

σ22[MPa]
58.40
46.69
34.98
23.27
11.56
-0.15

(a)λ=1.10(b)λ=1.20(c)λ=1.30

Gel phase 3

Fig.3.9 Change in distribution of tensile stress σ22 under tension (0◦ bunching
structure model).

32

3.3.330◦ 数珠繋ぎ構造モデルの解析結果
0.1

Experimentwith Gel phase 1
with Gel phase 2
with Gel phase 3

Hysteresis loss [J/m3]

Nominal Stress Σn22[MPa]

3

2

1

0.08

0.06

0.04

0.02

0

1

1.2

0

1.4

Stretch λ2(a) Nominal stress-stretch relations(b) Hysteresis loss

Fig.3.10 Comparison of (a) Nominal stress-stretch relations and (b) Hysteresis
loss of silica-filled rubber by simulation and experiment (30◦ bunching
structure model).

30◦ モデルにおいて界面ゲル相の物性を Gel phase 1 から Gel phase 3 まで 3 パターンに変えた解析から得られた (a) 公称応力―ストレッチ関係，(b) ヒステリシスロスを
実験データと比較して図 3.10 に示す．30◦ モデルの構造では，全体のストレッチ λ2 の増加に応じてシリカ粒子は左右に位置をずらすことが可能なため，ゲル相への応力集
中が 0◦ モデルに比べて少なくなる．そのためメッシュが計算途中で潰れることはなく，

1 サイクルの計算が可能であり，応力上昇挙動が実験結果に近づいた．ただし変形の集中するゲル相にアフィン変形を適用しているため，ヒステリシスロスは実験結果と比べてどの場合も小さくいずれのモデルもほぼ同じ値を示している．
最大ストレッチ λ2 = 1.50 における分子鎖ストレッチ λc の分布を図 3.11 に，引張方
向応力 σ22 の分布を図 3.12 に示す．30◦ モデルにおいても 0◦ モデルのときと同様，シリカ粒子間の距離が最小となる中心軸付近のゲル相に変形が集中し顕著な応力上昇を生じている．30◦ モデルの場合では粒子の側面のゴムマトリクス部の分子鎖ストレッチ λcは非対称となり，大きく延伸したゲル相側の変形が顕著となる．また (a) → (c) とゲル相の物性が軟らかいものになるにつれて粒子を繋ぎ止めている部分の分子鎖ストレッ33

チ λc は大きくなるが，図 3.4 で示したように変形抵抗が小さいため応力 σ22 は小さく
なる．

λc
2.75
2.40
2.05
1.70
1.35
1.00

(a)Gel phase1(b)Gel phase2(c)Gel phase3

Fig.3.11 Distribution of morecular chain stretch λc (30◦ bunching structure model).

σ22[MPa]
42.71
34.17
25.64
17.11
8.57
0.04

(a)Gel phase1(b)Gel phase2(c)Gel phase3

Fig.3.12 Distribution of tensile stress σ22 (30◦ bunching structure model).

34

3.3.445◦ 数珠繋ぎ構造モデルの解析結果
0.1

Experimentwith Gel phase 1
with Gel phase 2
with Gel phase 3

Hysteresis loss [J/m3]

Nominal Stress Σn22[MPa]

3

2

1

0.08

0.06

0.04

0.02

0

1

1.2

0

1.4

Stretch λ2(a) Nominal stress-stretch relations(b) Hysteresis loss

Fig.3.13 Comparison of (a) Nominal stress-stretch relations and (b) Hysteresis
loss of silica-filled rubber by simulation and experiment (45◦ bunching
structure model).

図 3.13 に 45◦ モデルの解析結果を示す．45◦ モデルの構造では，30◦ モデルの場合よりもさらに全体のストレッチ λ2 の増加に応じてシリカ粒子は左右に位置をずらすことが可能なため，系全体の応力上昇はさらに緩やかとなった．そのため実験結果とかなり近い応答を示している．ただしヒステリシスロスは実験結果と比べてどの場合も小さくほぼ同じ値を示しており，30◦ モデルの場合と同程度である．
最大引張時（λ2 = 1.50）における分子鎖ストレッチ λc の分布を図 3.14 に，引張方
向応力 σ22 の分布を図 3.15 に示す．これまでに説明した 0◦ モデル，30◦ モデルの解析結果と同様，45◦ モデルの場合でもシリカ粒子同士を繋ぎ止めているユニットセルの中心
軸付近のゲル相に変形が集中し顕著な応力上昇を生じている．また 30◦ モデルの場合と同様，粒子の片側側面部分のゴムマトリクス部が上下の粒子により引張変形を受けることで変形集中部のゲル相と同程度の分子鎖ストレッチ λc を生じている．この部分の応力 σ22 は 30◦ モデルに比べて高くなっているが，やはり粒子を連結している部分のゲル相の配向硬化が顕著であり系の応力上昇を担っているのは主としてゲル相である．

35

λc
2.36
2.09
1.81
1.54
1.27
1.00

(a)Gel phase1(b)Gel phase2(c)Gel phase3

Fig.3.14 Distribution of morecular chain stretch λc (45◦ bunching structure model).

σ22[MPa]
16.89
13.51
10.13
6.75
3.37
-0.01

(a)Gel phase1(b)Gel phase2(c)Gel phase3

Fig.3.15 Distribution of tensile stress σ22 (45◦ bunching structure model).

36ゲル相の非アフィンモデル導入

3.4.1

非アフィンモデルによるゲル相のモデル化

1.5

1

0.5

0

1

1.2

1.4

20wt.%
15wt.%
10wt.%
8wt.%
4wt.%
2wt.%
0wt.%

2

1.5

Nominal Stress σ22[MPa]20wt.%
15wt.%
10wt.%
8wt.%
4wt.%
2wt.%
0wt.%

2

Nominal Stress σ22[MPa]Nominal Stress σ22[MPa]

3.4

1

0.5

0

1

1.2

1.4

20wt.%
15wt.%
10wt.%
8wt.%
4wt.%
2wt.%
0wt.%

2

1.5

1

0.5

0

1

1.2

Stretch λ2

Stretch λ2

Stretch λ2(a)10%(b)15%(c)20%

1.4

Fig.3.16 Comparison of nominal stress-stretch relations of (a)10 ％, (b)15 ％,(c)20 ％ silica-filled rubber by experiment.

3.2.1 節で説明したゲル相物性見積もり実験の問題点として，シリカ充填ゴム中に実際に存在するゲル相を抽出して引張試験を行った結果ではなく，ゲル相の加硫度を評価し，それと同じ加硫度のゴムを引張試験した結果であることが挙げられる．図 3.16 にシリカ含有率を (a)10 ％，(b)15 ％，(c)20 ％ と変えて作成したシリカ充填ゴムの引張試験結果を示す．それぞれシリカ粒子とゴム材の界面状態を制御するカップリング剤の添加量を 0wt.％ から 20wt.％ まで変えて作成したシリカ充填ゴムで引張試験を行っており，安定した変形挙動が観察される 2 サイクル目の結果を比較している．カップリ
ング剤由来の硫黄はゴムマトリクス部では検出されておらず，界面ゲル相にのみその存在が実験的に確認されており，ゲル相の架橋密度はカップリング剤の量に依存する．

(a) のシリカ含有率 10 ％ の場合ではカップリング剤の量が多くなっても応力上昇挙動
およびヒステリシスにあまり変化はなく，またカップリング剤の添加量が 0wt.％ の場
合はどのシリカ含有率でも最大応力はほぼ同じでヒステリシスループの面積が僅かに増大している．一方シリカ含有率を増加させ，さらにカップリング剤の添加量を増大させた場合では顕著な応力上昇とヒステリシスの大幅な増大の傾向を示している．これは，シリカ含有率が低い場合ではシリカ粒子の周囲に生成するゲル相の量が少ない

37

Experiment
Ns=8.4 (Affine)Ns=8.4 (Non affine)Ns=6.4 (Non affine)Ns=4.4 (Non affine)

15

20

Segment N

True Stress σ22[MPa]20

Ns=8.4 (Affine)Ns=8.4 (Non affine)Ns=6.4 (Non affine)Ns=4.4 (Non affine)unloading

24

10

16

loading
12

unloading508loading

4
1

2

3

4

1

2

3

4

Stretch λStretch λ(a) True stress-stretch relations(b) Segment-stretch relations

Fig.3.17 Comparison of (a) True stress-stretch relations and (b) Segmentstretch relations by aﬃne model and non aﬃne model.

ため，カップリング剤によるゲル相の物性変化が系の応答に現れず，一方シリカ含有率とカップリング剤が共に増加すれば系全体に占めるゲル相の割合ならびに物性の寄与が大きくなるためと考えられる．3.2.2 節で説明したゲル相の変形過程にアフィン変形を導入したモデルではヒステリシスはほとんどないと仮定したが，図 3.16 はシリカ
含有率・カップリング剤が共に増加することでヒステリシスが大きくなることを示している．したがってここではゲル相のヒステリシス増加を表現するために，非アフィンモデルによるゲル相のモデル化を行う．

3.2.2 節で説明したアフィンモデルによる Gel phase 2 のゲル相と，8 鎖モデル A の要素 α に 2.3 節で説明した非アフィン分子鎖網目モデルを適用したゲル相について，

λ = 4.0 まで引張変形を与えた時の真応力ーストレッチ関係を図 3.17(a) に示す．Gelphase 2 の材料定数は全て同じまま非アフィン変形を適用した N s =8.4 のモデルに加えて，変形抵抗の増大についても検討するためゲル相の初期セグメント数 N s を N s =6.4，

4.4 に変化させたモデルを作成した．赤色の線が 3.2.1 節で説明した「ゲル相と同じ加硫度のゴム」の実験結果である．全てのゲル相モデルにおいて，総セグメント数

Nｎ (= nα Nα + nβ Nβ + nB NB ) は Nｎ = 1.95 × 1027 に統一している．各モデルのパラメータを表 3.3 に示す．N s =8.4（Aﬃne）モデルと N s =8.4（Non aﬃne）モデルを比較すると，初期セグメント数は同じであるが，非アフィンモデルは引張時にセグメン

38

ト数の増大（分子鎖同士の絡み点の減少）を生じて応力軟化する．そのためアフィン
モデルよりも応力上昇は少ないが，除荷時の応力低下が大きくなりヒステリシスは大きくなる．
図 3.17(b) に非アフィン変形による変形中のセグメント数変化を示す．アフィン変形の場合は変形中のセグメント数変化を生じないため負荷・除荷を通して N s =8.4 で一定値である．一方非アフィン変形を適用した場合はゲル相モデル全体の引張変形 λ2 = 4.0
に対するセグメント数の増加割合を 100 ％ としたため，初期セグメント数 N s =8.4，

6.4，4.4 に対して，λ2 = 4.0 において 2 倍の N s =16.8，12.8，8.8 となる．除荷時はセグメント数の変化は起こらない．
このような非アフィン変形を適用したゲル相の物性値を，3.3.1 節で検討したシリカ

2 粒子数珠繋ぎ構造で最も安定に解析を行うことのできた 45 °モデルに導入し，巨視座標系の x2 方向に u̇ = 100[mm/min] の一定な変形速度で最大伸びが λ2 = 1.5 になるまで変形を与えて 1 サイクルの変形挙動を解析する．その他の解析条件は 3.3 節と同様である．

Table 3.3 Parameter of rivised eight chain model (Gel phase, Non aﬃne model).

N s =8.4
CαR (= nα kB T )0.4283CβR = (nβ kB T )0.42

CBR (= nB kB T )
0.10

Nα
8.4

Nβ
8.4

NB
8.4

Ĉ1A
5.0 × 103

C2A
−0.50

mA
3.2

Ĉ1D
3.0 × 103

C2D
−0.50

mA
4.8

N s =6.4
CαR (= nα kB T )0.6447CβR = (nβ kB T )0.50

CBR (= nB kB T )
0.10

Nα
6.4

Nβ
6.4

NB
6.4

Ĉ1A
5.0 × 103

C2A
−0.50

mA
3.2

Ĉ1D
3.0 × 103

C2D
−0.50

mA
4.8

N s =4.4
CαR (= nα kB T )1.40

CβR = (nβ kB T )0.31

CBR (= nB kB T )
0.10

Nα
4.4

Nβ
4.4

NB
4.4

Ĉ1A
5.0 × 103

C2A
−0.50

mA
3.2

Ĉ1D
3.0 × 103

C2D
−0.50

mA
4.8

39

非アフィンゲル相モデルによる 45◦ 数珠繋ぎ構造の解析結果

3.4.2

1.5

Hysteresis loss [J/m3]

Nominal Stress Σn22[MPa]0.1Experiment
Ns=8.4 (Affine)Ns=8.4 (Non affine)

1

0.5

0

1

1.2

Stretch λ2

0.08

0.06

0.04

0.02

0

1.4

(a) Nominal stress-stretch relations(b) Hysteresis loss

Fig.3.18 Comparison of(a) Nominal stress-stretch relations and (b) Hysteresis
loss of silica-filled rubber by simulation and experiment (45◦ bunching
structure model).

初期セグメント数 N s =8.4 の非アフィンゲル相モデルを導入した 45◦ 数珠繋ぎ構造モデルの解析から得られた (a) 公称応力―ストレッチ関係，(b) ヒステリシスロスを前節での結果と比較して図 3.18 に示す．非アフィン変形を適用した場合では変形後期の応力上昇が抑制され最大引張時の応力は低下し，また実験結果に見られるような変形後期の硬化挙動（3 次曲線的な負荷曲線）は見られなくなった．これはゲル相内で非アフィン変形によるセグメント数 N s の増加（＝分子鎖同士の絡み点数の減少）を生じ，分子鎖の配向硬化が弱くなったためである．また除荷時の応力低下が顕著となりヒス
テリシスロスが増大しており，実験結果以上の数値を示している．
図 3.19〜図 3.21 に最大引張 (λ2 = 1.5) における分子鎖ストレッチ λc 分布，引張方
向応力 σ22 分布，セグメント数の変化量 ∆N 分布をそれぞれ示す．図 3.19 を見ると，ゲル相アフィン変形モデル (a) よりも，非アフィン変形モデル (b) の方が局所的に大きく延伸している分子鎖がシリカ粒子を連結している部分に認められる．図 3.21 のセグ
メント数の変化量 ∆N について，ゲル相アフィン変形モデルではゴム相に僅かに非アフィン変形を生じているだけであるが，ゲル相非アフィンモデルの図 (b) ではシリカ

40

粒子を連結する部分のゲル相に著しいセグメント数 N s の増加を生じている．本モデルにおいて，セグメント数 N s の増加は分子鎖のすり抜けによって絡み点数が減少することに対応しており，先述のように分子鎖の配向硬化を緩和し図 3.20(b) のように応力が低下する．
λc
2.64
2.31
1.98
1.66
1.33
1.00

(a)Affine(b)Non affine

Fig.3.19 Comparison of distributions of morecular chain stretch λc (45◦ bunching structure model).

σ22[MPa]
13.48
10.78
8.09
5.39
2.69
-0.01

(a)Affine(b)Non affine

Fig.3.20 Comparison of distributions of tensile stress σ22 (45◦ bunching structure model).

41

∆N
11.23
8.98
6.74
4.49
2.25
0.00

(a)Affine(b)Non affine

Fig.3.21 Comparison of distribution of change in the segment ∆N (45◦ bunching structure model).

ゲル相の初期セグメント数 N s を変化させた場合

3.4.3

0.15

Nominal Stress Σn22[MPa]2Hysteresis loss [J/m3]

Experiment
Ns=8.4Ns=6.4Ns=4.4

1.5

1

0.5

0

1

1.2

Stretch λ2

0.1

0.05

0

1.4

(a) Nominal stress-stretch relations(b) Hysteresis loss

Fig.3.22 Comparison of (a) Nominal stress-stretch relations and (b) Hysteresis
loss of silica-filled rubber by simulation and experiment (45◦ bunching
structure model, non aﬃne gel phase).

ゲル相の初期セグメント数 N s を変化させたときの (a) 公称応力―ストレッチ関係，

42

(b) ヒステリシスロスを図 3.22 に示す．λ2 = 1.5 まで延伸させた場合，N s =8.4 の場合に比べて N s =6.4，4.4 のようにセグメント数 N s が小さい（=分子鎖の絡み点数が多い）場合は，負荷時の分子鎖の配向硬化が強くなるため変形初期から応力が急激に上昇する．またそのような顕著な配向硬化によって負荷から除荷に移行した際の応力がより弾性的に低下しヒステリシスロスも増大する．したがって 3.4.1 節で説明したカッ
プリング剤の増加による変形抵抗とヒステリシスロスの増大の傾向をセグメント数 N sを減少させることにより再現可能であることが分かる．
図 3.23〜図 3.25 に最大引張 (λ2 = 1.5) 時の分子鎖ストレッチ λc 分布，引張方向応力

σ22 分布，セグメント数の変化量 ∆N 分布をそれぞれ示す．図 3.23 を見ると，(a) に比べて (b)，(c) とゲル相のセグメント数 N s が小さくなるほど粒子を繋ぎ止めている部
分の分子鎖ストレッチ λc は小さくなり，代わりに赤丸で示した粒子側面のゴムマトリ
クス部の分子鎖ストレッチ λc が大きくなっている．このため図 3.25 のセグメント数の変化量 ∆N 分布では (a) → (b) → (c) と初期セグメント数 N s が小さくなるほど変形集
中部での非アフィン変形は減少し，セグメント数の変化量 ∆N は小さくなる．図 3.24の応力分布を見ると，(a) → (b) → (c) と初期セグメント数 N s が小さくなるにつれゲ
ル相の変形抵抗が大きくなり，また図 3.25 に示したように変形集中部でのセグメント数の変化量 ∆N も小さくなるため，分子鎖の配向硬化が顕著になり変形集中部は非常に大きな応力 σ22 を示す．
λc
2.64
2.31
1.98
1.66
1.33
1.00

(a)N=8.4(b)N=6.4(c)N=4.4

Fig.3.23 Comparison of Distribution of morecular chain stretch λc (45◦ bunching structure model, non aﬃne gel phase).

43

σ22[MPa]
20.48
16.39
12.31
8.22
4.14
0.05

(a)N=8.4(b)N=6.4(c)N=4.4

Fig.3.24 Comparison of Distribution of tensile stress σ22 (45◦ bunching structure model, non aﬃne gel phase).

∆N
11.23
8.98
6.74
4.49
2.25
0.00

(a)N=8.4(b)N=6.4(c)N=4.4

Fig.3.25 Comparison of Distribution of change in the segment ∆N (45◦ bunching structure model, non aﬃne gel phase).

44

第4章
分子鎖セグメントの不均一分布が変形挙動に及ぼす影響

4.1

ゴム単相の変形挙動評価

4.1.1

解析モデル

Fig.4.1 Schematic view of heterogeneous distribution of molecular chain.

ゴム材に硫黄などの架橋剤を添加しゴム材を硫黄架橋した場合，図 4.1 に模式的に

45

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 (N0=14.0)

Segment Number N (N0=14.0)(a) Normal distribution(b) Double peaks distribution

Fig.4.2 Distribution pattern of segment numbers (Rubber matrix).

示すように，系内部に架橋点の多い（架橋点間のセグメント数 N が小さい）領域と少ない（セグメント数 N が大きい）領域が不均一に存在し，その大きさも場所によって大きく異なることが最近の研究により分かっている (5) ．このような硫黄架橋によるゴ
ム内部の分子鎖の絡み点の不均一性に関して，その分布形態を観察する研究 (6) が進められているが，実際のゴム内部における「セグメント数」の分布は不明である．そのため本章では非晶性高分子材料に存在する微視的な分子鎖密度の不均一性を表現した
シミュレーション (58) にならい，図 4.2 に示すように初期セグメント数 N を (a) 正規分
布と (b) ダブルピーク型分布としたときの応答について検討した．グラフの横軸は架橋点間の分子鎖セグメント数 N ，縦軸はユニットセル内の要素数 M である．ユニット
セル内の平均セグメント数 N0 = ΣN/ΣM を，3.1 節で説明したゴムマトリクス部のセグメント数 N と同じ N0 =14.0 に統一し，分散範囲を 0.75N0 から 1.25N0 としたもの
を ND1 および DP1，0.50N0 から 1.50N0 としたものを ND2 および DP2，0.25N0 から

1.75N0 としたものを ND3 および DP3 としている．
このようなセグメント数 N の不均一分布に従い，ゴム単相のユニットセルの要素に異なるセグメント数 N をランダムに割り当てたものを解析モデルとする．図 4.3 にダブルピーク型分布 DP3 に従う不均一性を導入したユニットセルのセグメント数 N の分布を例として示す．均質化法に基づく有限要素解析により，このような不均一分布をもつユニットセルが周期的に繰り返されるゴム単相の応答について，前章までの解46

析と同様の条件（平面ひずみ状態，u̇ = 100[mm/min]，最大伸びが λ2 = 2.0）で 1 サ
イクルの変形挙動を解析した．温度は変形過程を通して一定で T = 296[K] とし，ゴムの物性値は 3.1 節で説明したゴムマトリクス部のそれと同じにしている．

Ν
24.5
20.3
16.1

L0

x2

11.9
7.7

y2

x1

0

0

3.5

y1

Unit cell

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

4.1.2

ゴム単相モデルの解析結果

図 4.4，図 4.5 はそれぞれ正規分布，ダブルピーク型分布の不均一性を導入したときの(a) 公称応力―ストレッチ関係，(b) ヒステリシスロスである．図中の Homogeneousはセグメント数 N が全て N0 =14.0 の均一なゴムの計算結果（均質化有限要素法を使わずに構成式から計算される応答）である．分布を導入しない場合に比べ，引張中のゴムの応力上昇ならびに最大引張（λ2 = 2.0）時の応力は，正規分布，ダブルピーク型
分布のいずれもセグメント数 N の分散を大きくするほど大きな値を示した．ただし正規分布を導入した系では応力上昇はわずかであり，ダブルピーク型分布の方が顕著である．一方でヒステリシスロスはどの分布形態の場合でもほとんど変化がない．
図 4.6〜図 4.8 に各モデルでの最大引張時 (λ2 = 2.0) における分子鎖ストレッチ λcの分布，引張方向応力 σ22 分布，セグメント数の変化量 ∆N 分布をそれぞれ示す．図中の上部に記載した数字はユニットセル全体での各値（λc ，σ22 ，∆N ）の平均である．
図 4.6 の λc の分布について，ひずみ制御で引張を行っているので平均値 λ̄c はほぼ同じであるが，(a) → (c)，(d) → (f) と不均一分布の分散が大きくなるにつれ，分子鎖スト
レッチ λc の小さい部分（青〜青緑色）と分子鎖ストレッチ λc の大きい部分（黄緑〜赤色）のコントラストが明瞭となる．これはユニットセルに不均一分布を導入したこと

47

Homogeneous
ND1
ND2
ND3

0.8

0.12

Hysteresis loss [J/m3]

Nominal Stress Σn22[MPa]

1

0.6

0.4

0.2

0

1

1.5

0.09

0.06

0.03

0

2

Stretch λ2(a) Nominal stress-stretch relations(b) Hysterisis loss

Fig.4.4 Comparison of (a) nominal stress-stretch relations and (b) hysteresis loss of rubber (homogeneous vs. normal distribution of initial
segment number N ).
Homogeneous
DP1
DP2
DP3

0.8

0.12

Hysteresis loss [J/m3]

Nominal Stress Σn22[MPa]

1

0.6

0.4

0.06

0.03

0.2

0

0.09

1

1.5

0

2

Stretch λ2(a) Nominal stress-stretch relations(b) Hysterisis loss

Fig.4.5 Comparison of (a) Nominal stress-stretch relations and (b) Hysteresis
loss of rubber (homogeneous vs. double peaks distribution of initial
segment number N ).

で，絡み点の多い硬い部分（セグメント数 N の小さい部分）では変形が起こりにくく，反対に絡み点の少ない軟らかい部分（セグメント数 N の大きい部分）では大きな変形が生じるという，シリカ粒子を充填した時と同じ効果 (59) を生じるためである．ダブ
ルピーク分布を導入したモデルでは，正規分布を導入したモデルよりも硬い部分と軟らかい部分の差がより大きく，変形していない部分と大きく変形している部分がより48

明確に分かれている．
一方，図 4.7 の σ22 の分布は，応力―ストレッチ関係で述べたように不均一分布の分散が大きくなるにつれ，σ22 の平均値 σ̄22 は増加する．これも未充填ゴムに CB やシリカといった弾性体粒子を加えた時と同様の効果 (59) による．すなわち，変形しにくい
硬い部分（セグメント数 N の小さい部分）が増えると，系の中で変形を担うゴムの見かけの体積が減少し，分子鎖の配向硬化が強くなるため同じひずみ下では系全体の応力 σ22 が大きくなる．
ゴム単相の系では分子鎖ストレッチ λc に対するセグメント数の変化量 ∆N は非常に小さく，局所的な最大値で 1.0 にも達していない．そのため図 4.8 のセグメント数の変化量 ∆N の分布は，次節のゲル相の ∆N 分布のカラーバーの範囲を合わせたためすべて青となっており，∆N の不均一性はほとんど現れない．

49

λc=1.32

λc=1.32

λc=1.33(a)ND1(b)ND2(c)ND3λc=1.33

λc=1.33

λc=1.34(d)DP1(e)DP2(f)DP3

λc
2.28
2.04
1.81
1.57
1.33
1.09

Fig.4.6 Distribution of molecular chain stretch λc .

50

σ22=1.54σ22=1.57

σ22=1.64(b)ND1(c)ND2(d)ND3σ22=1.56

σ22=1.65

σ22=1.85(e)DP1(f)DP2(g)DP3

σ22[MPa]
5.86
4.99
4.12
3.25
2.38
1.51

Fig.4.7 Distribution of tensile stress σ22 .

51

∆N=0.15∆N=0.16∆N=0.16(a)ND1(b)ND2(c)ND3∆N=0.15∆N=0.16∆N=0.19(d)DP1(e)DP2(f)DP3

∆N
11.84
9.50
7.17
4.83
2.49
0.15

Fig.4.8 Distribution of change in the segment ∆N .

52

4.2

ゲル単相の変形挙動評価

4.2.1

解析モデル

前節のゴム単相の解析と同様に，系のセグメント数 N s に不均一分布を導入したゲル単相の有限要素解析を行う．図 4.9 にゲル単相の解析に用いたセグメント数 N s の分布を示す．ゲル相の物性値は 3.4.1 節で説明した非アフィン変形を導入した N s =8.4 の物
性値を用いるため各分布の平均セグメント数 N0s は N0s =8.4 としている．前節と同様，
分散範囲を 0.75N0s から 1.25N0s としたものを ND1 および DP1，0.50N0s から 1.50N0s としたものを ND2 および DP2，0.25N0s から 1.75N0s としたものを ND3 および DP3 と称
する．

500

Number of elements M

Number of elements M

500

400

300

200

100

400

300

200

100

0

0

0.25Ns0 0.50Ns0 0.75Ns0 Ns0 1.25Ns0 1.50Ns0 1.75Ns0
Segment number Ns (Ns0=8.4)

0.25Ns0 0.50Ns0 0.75Ns0 Ns0 1.25Ns0 1.50Ns0 1.75Ns0
Segment number Ns (Ns0=8.4)(a) Normal distribution(b) Double peaks distribution

Fig.4.9 Distribution pattern of segment numbers (Gel phase, N0s =8.4).

53

4.2.2ゲル単相モデルの解析結果
0.5

Homogeneous
ND1
ND2
ND3

1.5

Hysteresis loss [J/m3]

Nominal Stress Σn22[MPa]

2

1

0.5

0

1

1.5

0.4

0.3

0.2

0.1

0

2

Stretch λ2(a) Nominal stress-stretch relations(b) Hysteresis loss

Fig.4.10 Comparison of (a) Nominal stress-stretch relations and (b) Hysteresis
loss of gel phase (homogeneous vs. normal distribution of initial
segment number N ).
0.5

Homogeneous
DP1
DP2
DP3

1.5

Hysteresis loss [J/m3]

Nominal Stress Σn22[MPa]

2

1

0.5

0

1

1.5

0.4

0.3

0.2

0.1

0

2

Stretch λ2(a) Nominal stress-stretch relations(b) Hysteresis loss

Fig.4.11 Comparison of (a) Nominal stress-stretch relations and (b) Hysteresis
loss of gel phase (homogeneous vs. double peaks distribution of initial
segment number N ).

図 4.10，図 4.11 に得られた (a) 公称応力―ストレッチ関係，(b) ヒステリシスロスを示す．先のゴム相の場合とは異なり，小さな分散を与えた ND1 モデル，DP1 モデルでは最大引張時（λ2 = 2.0）の応力はごくわずかではあるが Homogeneous モデルよりも54小さくなった．その他のモデルではゴム単相の解析結果と同様，セグメント数 N s の
分散を大きくするほど引張時の応力上昇は大きくなり，ダブルピーク型分布で特に顕著である．またヒステリシスロスにも差を生じ，その差は大きくないものの，正規分布を導入したモデルでは ND3，ダブルピーク型の DP2 では微小に，DP3 モデルの場合は明白に増大している．
図 4.12〜図 4.14 に最大引張 (λ2 = 2.0) 時の分子鎖ストレッチ λc 分布，引張方向応力

σ22 分布，セグメント数の変化量 ∆N 分布をそれぞれ示す．図 4.12 の分子鎖ストレッチ λc 分布は，ゴム単相の場合と同様に，同じ平均値のもと不均一分布の分散が大きくなるにつれ分布のコントラストが強まる．なお，図 4.12(g) において最も変形が集中している部分の分子鎖ストレッチは λc = 2.28 となっており，ゴム単相解析の分子鎖スト
レッチ分布図 4.6(g) における変形集中部の分子鎖ストレッチ λc = 1.78 よりも大きい
数値となっている．同じ分子鎖ストレッチ λc に対するセグメント数の変化量 ∆N をゴ
ム相よりもゲル相の方が大きいと仮定しているため，局所的な ∆N の値がゴム単相の場合よりも非常に大きくなる（図 4.14）．このような変形集中部でのセグメント数 N sの増大（分子鎖の絡み点数の減少）は分子鎖の配向硬化を弱め変形抵抗が小さくなる．
したがって図 4.6(g) の変形集中領域では分子鎖の変形抵抗が減少し，周囲の硬い部分（セグメント数 N s の小さい部分）に引っ張られることでゴム単相解析の場合よりも局所的な分子鎖ストレッチ λc が大きくなる．
図 4.14 のセグメント数の変化量 ∆N の分布はゴム相のそれと比べ強いコントラス
トを示している．このような局所的なセグメント数の変化量 ∆N の増加が系全体の平均値 ∆N̄ の上昇をもたらし，ND3，DP2，DP3 モデルのヒステリシスロス増大につながったと考えられる．
なお，分布の分散の小さい ND1 や DP1 モデルでは平均値 ∆N̄ はゲル単相での値∆N =0.95 より大きくなっているにもかかわらず，先述のように僅かに最大引張時の応力は低下しておりヒステリシスロスも小さくなっている．これらのモデルでは，変形後期における局所的なセグメント数 N s の増加による応力緩和が，硬い部分の割合による見かけの体積減少による応力上昇を上回ったためであると考えられる．

55

λc=1.32

λc=1.33

λc=1.33(a)ND1(b)ND2(c)ND3λc=1.32

λc=1.33

λc=1.35(d)DP1(e)DP2(f)DP3

λc
2.28
2.04
1.81
1.57
1.33
1.09

Fig.4.12 Distribution of molecular chain stretch λc .

56

σ22=2.86

σ22=2.89

σ22=3.00(a)ND1(b)ND2(c)ND3σ22=2.86

σ22=3.02

σ22=3.35(d)DP1(e)DP2(f)DP3

σ22[MPa]
5.86
4.99
4.12
3.25
2.38
1.51

Fig.4.13 Distribution of tensile stress σ22 .

57

∆N=1.27∆N=1.30∆N=1.37(a)ND1(b)ND2(c)ND3∆N=1.29∆N=1.37∆N=1.62(d)DP1(e)DP2(f)DP3

∆N
11.84
9.50
7.17
4.83
2.49
0.15

Fig.4.14 Distribution of change in the segment ∆N .

58

4.3

ゲル単相の平均セグメント数 N0s を変化させた場合

4.3.1

解析モデル

前節と同様のセグメント数 N s の不均一分布を導入したゲル単相の系について，平均セグメント数 N0s を N0s =6.4，4.4 と変化させた解析を行う．ゲル相の物性値には，3.4.3
節の非アフィンモデルの N s =6.4，4.4 の物性値を用いる．前節で対象とした分散形態

ND1〜ND3，DP1〜DP3 について解析した．

4.3.2

N0s =6.4,4.4 としたゲル単相モデルの解析結果

図 4.15〜図 4.18 に得られた結果をまとめて示す．各図で公称応力ーストレッチ関係のグラフの縦軸のスケールが 4.2.2 節の解析結果の図と異なっていることに注意されたい．なお，ヒステリシスロスの図のスケールは 4.2.2 節と同じにしている．λ2 = 2.0 までの引張変形では，平均セグメント数 N0s を小さくするほど Homogeneous モデルの最大応力・ヒステリシスロスが増加しており，さらに前節で説明した不均一分布の分散を大きくした場合の応力・ヒステリシスロスの増加傾向がより顕著になっている．
図 4.19 に不均一分布による変化が最も大きい DP3 モデルのセグメント数の変化量∆N の分布を示す．上に示した ∆N̄ は平均値であり，均一な系での計算では N0s =8.4，
6.4，4.4 それぞれで ∆N =0.95，0.72，0.49 と小さくなっている．∆N が小さいほど除荷時の応力低下は大きくなるが，図 4.10，図 4.15，図 4.17 の均一な系のヒステリシス
ロスを見ると N0s が小さくなるほど大きくなっている．これは 3.4.4 節で述べたように，
ゲル相のセグメント数 N0s を小さくすると配向硬化が顕著になるため負荷時の応力上昇・除荷時の応力低下の勾配が急になり（＝応力ーストレッチ曲線の「長さ」が大きくなり），少ないセグメント数の変化量 ∆N も大きなヒステリシスロスを生じるようになるためである（ セグメント数の変化量 ∆N に対するヒステリシスロスの発現量が大きくなる）．一方，不均一分布を導入した系での計算では局所的に ∆N が増大する
ため，図 4.19 において均一な系に比べて ∆N̄ は大きくなっており，除荷時の応力低下
sが大きくなりヒステリシスロスが増大する．したがって，⃝
1 初期セグメント数 N0 が小さいほど応力上昇によるヒステリシスロスの発現量が大きくなり，かつ⃝
2 不均一分布を導入すると ∆N が増大するため，この両者の相乗効果によってヒステリシスロスは増大する．

59

0.5

Homogeneous
ND1
ND2
ND3

4

Hysteresis loss [J/m3]

Nominal Stress Σn22[MPa]

5

3

2

1

0

1

1.5

0.4

0.3

0.2

0.1

0

2

Stretch λ2(a) Nominal stress-stretch relations(b) Hysteresis loss

Fig.4.15 Comparison of (a) Nominal stress-stretch relations and (b) Hysteresis
loss of gel phase (N0s =6.4, homogeneous vs. normal distribution of
initial segment number N ).

0.5

Homogeneous
DP1
DP2
DP3

4

Hysteresis loss [J/m3]

Nominal Stress Σn22[MPa]

5

3

2

1

0

1

1.5

0.4

0.3

0.2

0.1

0

2

Stretch λ2(a) Nominal stress-stretch relations(b) Hysteresis loss

Fig.4.16 Comparison of (a) Nominal stress-stretch relations and (b) Hysteresis
loss of gel phase (N0s =6.4,homogeneous vs. double peaks distributionof initial segment number N ).

60

0.5

Homogeneous
ND1
ND2
ND3

4

Hysteresis loss [J/m3]

Nominal Stress Σn22[MPa]

5

3

2

1

0

0.4

0.3

0.2

0.1

1

1.5

0

2

Stretch λ2(a) Nominal stress-stretch relations(b) Hysteresis loss

Fig.4.17 Comparison of (a) Nominal stress-stretch relations and (b) Hysteresis
loss of gel phase (N0s =4.4, homogeneous vs. normal distribution of
initial segment number N ).

0.5

Homogeneous
DP1
DP2
DP3

4

0.4

Hysteresis loss [J/m3]

Nominal Stress Σn22[MPa]

5

3

2

1

0

1

1.5

0.3

0.2

0.1

0

2

Stretch λ2(a) Nominal stress-stretch relations(b) Hysteresis loss

Fig.4.18 Comparison of (a) Nominal stress-stretch relations and (b) Hysteresis
loss of gel phase (N0s =4.4, homogeneous vs. double peaks distributionof initial segment number N ).

61∆N=1.62∆N=1.28∆N=0.92(a)Ns=8.4(b)Ns=6.4(c)Ns=4.4

∆N
11.84
9.50
7.17
4.83
2.49
0.15

Fig.4.19 Distribution of change in the segment ∆N (DP3 model).

62

第5章
分子鎖の不均一性を反映させた構成
式の構築
5.1

セグメント不均一性を反映したゴム単相・ゲル単相の構成式化

Nominal Stress Σn22[MPa]245c

An element
Unit cell model
1.5

1

0.5

0

1

1.5

2

Stretch λ2(a)Silica filled-rubber(b)Silica filled-rubber model(c)Heterogeneous model(d)Nominal stress-stretch relations

Fig.5.1 Hierarchical model of silica-filled rubber.

第 3 章では，シリカ充填ゴムの局所的な構造（図 5.1(a)）を単純化して図 5.1(b) に示す直鎖状数珠繋ぎ構造モデルの解析を行った．第 4 章ではシリカ充填ゴム中のゴム
マトリクス部およびゲル相について，実験によって示唆される分子鎖の不均一分布を想定したモデル（図 5.1(c)）を対象に様々な検討を行った．図 5.1(d) は，分子鎖セグメ
ントの不均一分布を導入したゴム単相モデルの解析から得られた公称応力ーストレッ
チ関係であり，灰色で示した多数の曲線はユニットセル中の各要素の応答，赤色がこれらの要素からなる周期セルの平均化された巨視的な応答である．本節では，このようなゴム相・ゲル相の不均一性による応答を図 5.1(b) のようなシリカ充填構造の応答63評価に導入することを検討する．しかしながら，分子鎖の不均一性とシリカ粒子のス
ケールには大きなギャップがあり，2 段階のスケールアップを均質化法によって直接扱うことは現時点では困難である．そこで，ゴム単相・ゲル単相について行った前章の解析結果（図 5.1(d) の赤線）を，図 5.1(c) の各要素メッシュの計算に用いたものと同様の構成式モデル（図 5.2）でスケールアップしてゴム相・ゲル相のパラメータを決定
し，図 5.1(b) のシリカ充填構造の粘弾性応答を議論する．

BA
Ș

Affine8-chain modelșAffine model

Nonaffine
model

D

Ț

Viscous damper
Viscous damper

Fig.5.2 Revised eight chain model (nonaﬃne).

ゴム単相の構成式化

Nominal Stress Σn22[MPa]5.1.1

Heterogeneous model (DP3)
Revised 8 chain model

1.5

1

0.5

0

1

1.5

2

Stretch λ2

Fig.5.3 Nominal stress-stretch relations of rubber matrix by simulation (100[mm/min]).

64

Heterogeneous model (DP3)
Revised 8 chain model

0.3

Hysteresis loss[J/m3]

Maximam nominal stress Σn22 max[MPa] at λ2=2.0

Heterogeneous model (DP3)
Revised 8 chain model

2.4
2.1
1.8
1.5
1.2

0.2

0.1

0.9
10-3

10-2

10-1

100

101

102

0
10-3

103

Strain rate ε [1/s](a)Max nominal stress-strain rate relations

10-2

10-1

100

101

.
Strain rate ε [1/s]

102

103

(b)Hysteresis loss-strain rate relations

Fig.5.4 (a)Max nominal stress-strain rate relations and (b) hysteresis lossstrain rate relations of rubber matrix by simulation changed the strain
rate ε̇.

4.1 節で説明した，各メッシュのセグメント数 N の不均一分布を直接扱ったゴム単相の解析結果のうち，不均一分布の導入の影響が大きかった DP3 の結果を基に，図 5.2の非アフィン粘弾性 8 鎖モデルのパラメータフィッティングを行った．得られた構成
式のパラメータを表 5.1 に示す．3.1 節のゴムマトリクス部と同様に，総セグメント数

Nｎ (= nα Nα + nβ Nβ + nB NB ) は Nｎ = 1.95 × 1027 としている．解析条件は平面ひず
み状態，変形速度 u̇ = 100[mm/min]，最大伸びが λ2 = 2.0，温度 T = 296[K] である．
フィッティングにより得られた公称応力ーストレッチ関係を図 5.3 に示す．得られた「マクロな」構成式と，均質化法で不均一性を実際に考慮したシミュレーションそれぞれで，ひずみ速度を ε̇ = 1.0 × 10−3 〜1.0 × 103 [1/s] と変化させて最大伸び λ2 = 2.0 ま
での繰り返し変形を行った時の最大応力とヒステリシスロスを，横軸にひずみ速度をとって図 5.4(a)(b) に示した．最大応力のひずみ速度依存性は一致しており，ヒステリ
シスロスについてはずれを生じるものの，ヒステリシスロスが最大となるひずみ速度が ε̇ = 1.0 × 101 [1/s] で一致し，全体的な傾向も良好に再現している．

65

Table 5.1 Parameter of rivised eight chain model (Rubber matrix based on the
result of heterogeneous model (DP3)).

CαR (= nα kB T )0.2666CβR = (nβ kB T )0.3043

CBR (= nB kB T )0.1217

Nα
11.5

Nβ
11.5

NB
11.5

Ĉ1A
5.0 × 105

C2A
−0.50

mA
3.2

Ĉ1D
3.0 × 105

C2D
−0.50

mA
4.8

5.1.2

ゲル単相の構成式化

前節と同様に，前章のゲル相の DP3（初期セグメント数は N s =8.4）の結果を基にパラメータフィッティングを行った．公称応力ーストレッチ関係を図 5.5 に，得られたパラメータを表 5.2 に示す．総セグメント数 Nｎ は同じく Nｎ = 1.95 × 1027 としている．またひずみ速度を ε̇ = 1.0 × 10−3 〜1.0 × 103 [1/s] まで変化させて λ2 = 2.0 までの変形シミュレーションを行い，それぞれの (a) 最大応力とひずみ速度の関係，(b) ヒス
テリシスロスとひずみ速度の関係を図 5.4 に示した．最大応力のひずみ速度依存性は一致しており，(b) はやはり実際にセグメント数 N s の不均一性を考慮した計算からはずれを生じるものの，全体的な傾向は一致しておりヒステリシスが最大となるひずみ

Nominal Stress Σn22[MPa]速度も ε̇ = 1.0 × 100 [1/s] で一致している．

Heterogeneous model (DP3)
Revised 8 chain model

2.5

2

1.5

1

0.5

0

1

1.5

2

Stretch λ2

Fig.5.5 Nominal stress-stretch relations of gel phase by simulation (100 [mm/min]).

66

Heterogeneous model (DP3)
Revised 8 chain model

0.3

Hysteresis loss[J/m3]

Maximam nominal stress Σn22 max[MPa] at λ2=2.0

Heterogeneous model (DP3)
Revised 8 chain model

2.4
2.1
1.8
1.5
1.2

0.2

0.1

0.9
10-3

10-2

10-1

100

101

102

0
10-3

103

10-2

10-1

100

101

102

103

Strain rate ε [1/s]Strain rate ε [1/s](a)Max nominal stress-strain rate relations(b)Hysteresis loss-strain rate relations

Fig.5.6 (a)Max nominal stress-strain rate relations and (b) hysteresis lossstrain rate relations of gel phase by simulation changed the strain
rate ε̇.

Table 5.2 Parameter of rivised eight chain model (Gel phase based on the resultof heterogeneous model (DP3)).

CαR (= nα kB T )0.5214CβR = (nβ kB T )0.5113

CBR (= nB kB T )0.1217

Nα
6.9

Nβ
6.9

NB
6.9

Ĉ1A
5.0 × 103

C2A
−0.50

mA
3.2

Ĉ1D
3.0 × 103

C2D
−0.50

mA
4.8

67

不均一性を反映した構成式による 45◦数珠繋ぎ構造モ

5.2

デルの解析

0.15

Experiment
Homogeneous model
Pseudo-heterogeneous model

1.5

Hysteresis loss [J/m3]

Nominal Stress Σn22[MPa]

2

1

0.5

0

1

1.2

0.1

0.05

0

1.4

Stretch λ2(a) Nominal stress-stretch relations(b) Hysteresis loss

Fig.5.7 Comparison of(a) Nominal stress-stretch relations and (b) Hysteresis
loss of silica-filled rubber by simulation and experiment (45◦ bunching
structure model).

前節で作成したゴムマトリクス部およびゲル相の構成式を用いて 45◦ 数珠繋ぎ構造
モデルの変形シミュレーションを行った時の (a) 公称応力―ストレッチ関係，(b) ヒス
テリシスロスを図 5.7 に示す．図中黒線は 3 章の結果で，赤線は実験結果である．不均一性を導入するとゴム相・ゲル相どちらも最大応力・ヒステリシスロスともに増加するのでいずれも 3 章のそれより増大している．このため実験結果からさらに外れた結果となっているが，45◦ 数珠繋ぎ構造と実際のシリカ粒子のランダム分布の差もあるた
めここでは問題としない．
図 5.8〜図 5.10 に最大引張 (λ2 = 1.5) における分子鎖ストレッチ λc 分布，引張方向
応力 σ22 分布，セグメント数の変化量 ∆N 分布をそれぞれ示す．図では 3 章の結果と比較して示している．図中上部に記載している数字はユニットセル内のゴムマトリクス部 (R)，ゲル相 (G) それぞれにおける λc ，σ22 ，∆N の平均である．図 5.8 の分子鎖ス
トレッチ λc 分布は不均一性を考慮した場合もほとんど変化がないが，図 5.9 の引張方

68

向応力 σ22 分布では，粒子を繋ぎ止めているゲル相の変形集中領域とその周囲のゴムマトリクス部で分子鎖の配向硬化が顕著になり応力が上昇している．一方，図 5.10 のセグメント数の変化量 ∆N 分布では，ゴム部の ∆N̄ R は変わらないがゲル相の平均値∆N̄ G は僅かに増加している．これは，4.2.2 節で述べたゲル相中の軟らかい部分（セグメント数 N の大きい部分）での局所的な非アフィン変形によってゲル相全体の平均
セグメント数変化量 ∆N̄ が増加しているためである．このような変形集中部のゲル相における顕著な応力 σ22 の上昇とセグメント数の変化量 ∆N の増大が，図 5.7(b) に示
すヒステリシスロスの増大をもたらしたと考えられる．

λc

λRc=1.23

λRc=1.23

λGc=1.49

λGc=1.48

2.65
2.32
1.99
1.66
1.33
1.00

(a)Homogeneous model(b)Pseudo-heterogeneous model

Fig.5.8 Comparison of distributions of morecular chain stretch λc (45◦ bunching structure model).

69

σ22[MPa]σR22=1.38

σR22=1.58

σG22=3.54σG22=3.89(a)Homogeneous model(b)Pseudo-heterogeneous model

8.15
6.53
4.91
3.28
1.66
0.04

Fig.5.9 Comparison of distributions of tensile stress σ22 (45◦ bunching structure model).

∆N∆NR=0.12∆NR=0.11∆NG=2.59∆NG=2.78

13.11
10.49
7.86
5.24
2.62
0.00

(a)Homogeneous model(b)Pseudo-heterogeneous model

Fig.5.10 Comparison of distribution of change in the segment ∆N (45◦ bunching structure model).

70

第6章
結論
シリカ充填ゴムは，CB 充填ゴムとは異なる粒子分布形態やゴム部と界面との相互作用により幅広い力学特性を発現する可能性を秘めている．本研究では，実験により確認されたシリカ充填ゴムの数珠繋ぎ構造についての基礎的な検討として，種々の実験事実を反映させたシミュレーションモデルの構築，ならびに充填構造やゲル相の物性が巨視的な力学的特性に与える影響を検討した．以下に得られた結果をまとめて示す．
第 2 章では分子鎖網目理論に基づく粘弾性 8 鎖モデルを概説し，ゴム粘弾性体の構成式について説明した．さらに，粘弾性構成式を絡み点数が変形とともに変化することを許容する非アフィン分子鎖網目理論により一般化し，更新ラグランジュ法に基づく有限要素均質化法を定式化した．
第 3 章では，数珠繋ぎ構造を呈するシリカ充填ゴムの局所的な構造を模擬したユニットセルモデルを作成し，粒子の配置形態，ならびにゲル相の特性変化がシリカ充填ゴムの応答に及ぼす影響を検討した．ゴム相，ゲル相の実験に基づきモデル化した構成式を用いて，ユニットセル中の 2 個のシリカ粒子の配置を 0◦ ，30◦ ，45◦ に変化させた系で変形シミュレーションを行った．ただし，ゲル相にはセグメント数変化を生じないアフィンモデルを適用している．ユニットセル中に直列（0◦ ）にシリカ粒子を配置した場合，ユニットセルの中心軸付近の粒子を繋ぎ止めているゲル相に変形が集中して計算不能となったが，シリカ粒子配置が 30◦ ，45◦ の系では，系全体の変形 λ2 に応じてシリカ粒子は位置を左右にずらすことが可能であるため応力上昇は緩やかになり，1サイクルの応答を得ることが可能であった．変形後期に強い配向硬化を示すアフィン
ゲル相モデルの導入により，実験結果にみられる変形後期の顕著な応力上昇（S 字状の応力ーストレッチ曲線）も再現可能であることを示した．ただし，アフィンモデル
を仮定した解析ではヒステリシスロスに変化がない．そこで，非アフィンモデルによ71

るゲル相のモデル化を行い，数珠繋ぎ構造の応答に与える影響を評価した．初期セグ
メント数 N s の同じアフィンゲル相での解析結果と比較すると，変形後期においてゲ
ル相内のセグメント数 N s が増加し，分子鎖の配向硬化が弱くなるため系の最大応力は低下するが，除荷時の応力低下も大きくなりヒステリシスロスが増大することが示された．さらに，ゲル相の初期セグメント数 N s の効果について検討し，初期セグメ
ント数 N s が小さいほど負荷時の応力上昇挙動ならびに負荷から除荷に移行した際の応力低下がより顕著になり，ヒステリシスロスの増大をもたらすことを示した．これ
により，非アフィンゲル相モデルを導入した数珠繋ぎ構造モデルにおいて，N s を小さくすることにより実験結果に見られるカップリング剤の増加に伴う変形抵抗とヒステリシスロスの増大の傾向を再現可能であることを示した．
第 4 章では，実験により報告されているゴムを硫黄架橋した際に生じる絡み点の不均一性について，ゴム相・ゲル相の各有限要素メッシュにセグメント数 N をランダムに変化させたモデルにより検討した．いずれも，不均一分布を導入することで，変形しにくい硬い部分（セグメント数 N の小さい部分）が，系の中で変形を担うゴム・ゲ
ルの見かけの体積減少（＝ゴムに CB やシリカ等の弾性体粒子を加えた時と同様の効果）をもたらすため，最大応力が上昇する．ゲル相ではさらに，変形しやすい軟らかい部分（セグメント数 N s の大きい部分）への変形集中・非アフィン変形の発生により，局所的な ∆N の値が増大し，除荷時の応力を低下させヒステリシスロスが増大する．さらに，ゲル単相の平均セグメント数 N0s を小さくする（硬いゲル相を想定する）と，3 章で述べた応力上昇・ヒステリシスロスの増大が，局所的な不均一性によって強められることを示した．
第 5 章では，第 4 章の解析から得られた絡み点の不均一性をもつゴム相・ゲル相の応答を，粘弾性 8 鎖モデルによる構成式で再度パラメータフィッティングしてスケー
ルアップを行い，第 3 章の 45◦ 数珠繋ぎ構造の解析を行って，絡み点の不均一性がシリ
カ充填ゴムの応答に及ぼす影響を評価した．その結果，不均一性の導入によりゴム相・ゲル相共に分子鎖の配向硬化が顕著になり，第 3 章で解析した不均一性をもたないモデルの結果よりも最大応力・ヒステリシスロスともに大きくなることを示した．
本研究で用いたシリカ充填ゴムの直鎖状数珠繋ぎ構造モデルは，実際のシリカ充填
ゴムの局所的な一部分を取り出し単純化を行ったモデルであるため，実際のシリカ充
填ゴムの自由度を大きく制限したものとなっており，実験結果に見られる以上に顕著な応力上昇および大きなヒステリシスロスを示している．実験結果における変形後期

72

の顕著な応力上昇挙動は，シリカ粒子の凝集構造の変化に起因する可能性も否定できない．今後，実際のシリカ充填ゴムに見られる巨視的初期等方性応答を呈する，より大きなスケールのモデルによる検討が待たれる．

73

参考文献
(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) 柴山充弘, 日本ゴム協会誌, 84-1,(2011),14-20.
(6) 中島健, 岩蕗仁, 浦部匡史, 伊藤万喜子, 藤波想, 西敏夫, 日本ゴム協会誌, 857,(2012),46.
(7) 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.
(8) 池田裕子, シンクロトロン X 線による架橋天然ゴムの伸長結晶化解析, 日本ゴム協会誌, 84-1,(2011),pp.29-36(9) Gao, J. and Weiner, J. H., Computer Simulation of Viscoelasticity in Polymer
Melts, Macromolecules, 25, (1992), 1348-1356.
(10) 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.
(11) Ying, Li, Martin, Kroger. and Wing, Kam, Liu., Primitive chain network study
on uncrosslinked and crosslinked cis-polyisoprene polymers, Polymer, 52, (2011),
5867-5878.
74

(12) 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 andScience, CDROM(13) 屋代如月, 内藤正登, 皆川康久, 冨田佳宏, 非晶性高分子材料のヒステリシスに関する分子動力学的研究, 日本機械学会論文集, A72, (2006), 277-284.

(14) 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.
(15) . 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.
(16) 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.
(17) 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.
(18) 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.
(19) 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.
(20) F. Bueche., Molecular Basis for the Mullins Eﬀect, Journal of Applied Polymer
Science, Vol.4, (1960), pp.107-114.

75

(21) 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.
(22) 内藤陽子, 伊藤眞義, シリカ充てん加硫ゴム材料の応力軟化に関する研究, 日本ゴム協会誌, 82-9,(2009),pp.394-399(23) 古谷泰大, 陸偉, 内藤正登, 冨田佳宏, 分子鎖網目モデルによるカーボンブラック
充填ゴムの粘弾性変形挙動シミュレーション, 日本機械学会第 17 回計算力学講演会講演論文集, No.04-40(2004-11), pp.17-18.

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

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

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

76

(33) Kuhn, W. and Grun, F., Beeziehuugen Zwischen Elastischen Konstanten und
Dehuungsdoppelbrechung Hochelastischer Stoﬀe, Kollooidzeitschrift, Vol.101,(1942), pp.248-271.
(34) 土井正男, 小貫明, 高分子物理・相転移ダイナミクス, (1992), 岩波書店.
(35) James, H.M. and Guth, E., Theory of the Elastic Properties of Rubber,J.Chem.Phys., 11 (1934), 455-481.
(36) James, H.M. and Guth, E., Theory of the elastic properties of rubber, Journal ofChemical Physics, Vol.11, (1943), pp.455-481.
(37) 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.
(38) 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.
(39) 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.
(40) 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.
(41) Green, A.E. and Zerna, W., Theoretical Elasticity, (1968), Oxford University
Press.
(42) Tomita,Y.,Azuma,K.and Naito,M.,Int.J.Mech.Sci., Vol50, No.5,(2008),pp856868.
(43) Truesdell, C. and Noll, W., The Non-Linear Field Theories of Mechanics,
Springer-Verlag, (1965).
(44) 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.
(45) 北川浩, 弾・塑性力学−非線形解析のための基礎理論−, (1987), 裳華房.
77

(46) 久田俊明, 非線形有限要素法のためのテンソル解析の基礎, (1992), 丸善.
(47) 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.
(48) Doi, M. and Edwards, S.F., The Theory of Polymer Dynamics, (1986), Oxford
University Press.
(49) Gennes, P.G. de., Reptation of a Polymer Chain in the Presence of Fixed Obstacles, Journal of Chemical Physics, Vol.55, (1971), pp.572.
(50) Raha, S. and Bowden, P.B., Birefringence of Plastically Deformed Polymethyl
Methacrylate, Polymer, Vol.13, (1972)，pp.174-183.
(51) 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.
(52) 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.
(53) 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.
(54) 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.
(55) 比嘉吉一, 冨田佳宏, 粒子強化型複合材の均質化法による変形挙動のモデル化とシ
ミュレーション, 日本機械学会論文集Ａ編, Vol.66, No.648(2000), pp.1441-1446.

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

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

11), pp.397-398.

79

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

1
σi =
3

{

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

式 (A.1) を左 Cauchy-Green 変形テンソル Aij を用いて表すと，
(
)}
(
{ √
)
√
λc
1 −1
1
λβ
λ2i
−1
R
R
√
√
σij =
Cα Nα L+Cβ Nβ L
3
λγ
λc
Nα
Nβ
{
(
)}
′
√
λcB
λ2i
1
R
+
CαB
NαB L−1 √− 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

80

(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β λ̇γ
λβ ζ ′
1
ζ
L′′
′
R
L +√
Aij +
CαB NαB √−A′ij A′kl /A′mm− 2√
3
λγ 3AmmNαB λcBNβ√]R}
L′′ CαBNαB {+δik A′jl + A′ik δjl (ε̇kl − ε̇pkl ) − ṗδij
(A.4)
λcB
∇

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

81

√
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
82

(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 y11
2Ae1 = det  11
1
2Ae3 = det 1
1

,

,

x
x2
x3
x
x1
x2


y
y2 
y3

y
y1 
y2

(B.5)

(B.6)

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

ψ1 =

Ae1
Ae

,

ψ2 =

Ae2
Ae

,

ψ3 =

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

83

Ae3
Ae

(B.7)

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

付録 C
関連発表論文・講演論文
学術講演 中田伸哉 ，屋代如月，冨田佳宏，FEM シミュレーションによるシリカ充填ゴムの履歴現象の評価，日本機械学会 第 24 回計算力学講演会 CMD2011, 岡山大学,(2011.10)
 Y. Tomita, S. Nakata, M. naito, K. Yashiro, Evaluation of Deformation Induced
Texture Dependent Mechanical Characteristics of Silica-Filled Rubber, International Symposium on Plasticity 2012 and Its Current Applications, San Juan,
Puerto Rico, (2012.1)
 Y. Tomita, S. Nakata, M. naito, K. Yashiro, Evaluation of Micro- to Macroscopic Deformation Behavior of Silica-Filled Rubber under Monotonic and Cyclic
Straining, WCCM2012, The 10th World Congress on Computational Mechanics,
Sao Paulo, Brazil, (2012.7)
 中田伸哉，屋代如月，冨田佳宏，FEM シミュレーションによるシリカ充填ゴムの微視構造のモデル化と力学特性評価，日本機械学会 第 25 回計算力学講演会

CMD2012, 神戸ポートアイランド南地区, (2012.10)
 Y. Tomita, S. Nakata, M. naito, K. Yashiro, Evaluation of Deformation Behavior of Silica-Filled Rubber under Monotonic and Cyclic Straining, AEPA2012,
Asia-Pacific Conference on Engineering Plasticity and its Applications, Nanyang
Technological University, Singapore, (2012.7)84論文貼り付け用ページ85

論文貼り付け用ページ86論文貼り付け用ページ87

論文貼り付け用ページ

88論文貼り付け用ページ

89

論文貼り付け用ページ90論文貼り付け用ページ

91論文貼り付け用ページ92論文貼り付け用ページ

93論文貼り付け用ページ

94論文貼り付け用ページ

95論文貼り付け用ページ96

謝辞
本研究を遂行するにあたり，著者が大学在籍時より研究全般に渡り懇切丁寧な御教授，御指導を賜るとともに，貴重な研究発表の機会や素晴らしい研究環境を提供して頂いた屋代如月准教授に深甚な感謝の意を表します．また，日頃から研究について多くの御助言，御教示を賜るとともに，本論文の作成にあたり丁寧かつ熱心なご指導を賜りました冨田佳宏教授（現福井工業大学）に心より感謝致します．
本研究に対して貴重な御助言と御指導，有益な議論を賜りました田中克志教授に深く感謝致します．また，研究活動が円滑に行えるよう数々の便宜をはかって頂きました古宇田由夫技術職員に心より感謝の意を表します．
本研究を進めるうえで，貴重な御助言と数々の便宜をはかって頂きました内藤正登
氏（住友ゴム工業株式会社）をはじめ，本研究において実験データを提供していただ
いた住友ゴム工業株式会社の諸氏に厚く御礼申し上げます．
同じ研究テーマに取り組み，有益な討論と適切なご指摘を頂いた，望月利紀氏（現ヤマハ発動機株式会社）に心から感謝致します．また，後輩として活発な議論を交わした本馬剛徳氏に感謝致します．
日々の研究活動や，大学院 1 年次の就職活動中に多くの有益なご意見を頂いた表面・界面工学研究室の諸先輩方に深く感謝いたします．また研究室配属以来，同輩として互いに切磋琢磨し，時には励ましあった池宮一繁氏，住谷昂大氏，竹下和也氏，西村
英晃氏，藤原正大氏に厚く御礼申し上げます．その他の表面・界面工学研究室諸氏ならびに，固体力学研究室諸氏にも今日に至るまで数々の御支援と御協力を頂きましたので，ここに記し感謝の意を表します．
最後に，6 年間の長きに渡る学生生活において，多大な支援と理解を頂いた両親に，敬意と感謝の意を表します．
平成 25 年 2 月
中田97

伸哉