(57)【要約】    【課題】  ゴム材料の実際の変形挙動に近づけたシミュレーションを行う。    【解決手段】  マトリックスゴムとフィラーとを含むゴム材料の変形をシミュレーションするゴム材料のシミュレーション方法であって、少なくとも一つの前記フィラーを有限個の要素に分割したフィラーモデルと、前記マトリックスゴムを、初期状態で前記フィラーモデルの外周面と一致ししかも前記外周面から離間可能な内周面を有した有限個の要素に分割したマトリックスモデルとを含むゴム材料モデルを設定するステップ、前記ゴム材料モデルに引張変形の条件を与えて変形計算を行うステップ、及び前記ゴム材料モデルから必要な物理量を取得するステップを含み、かつ、前記変形計算を行うステップは、前記マトリックスモデルの内周面と、前記フィラーモデルの外周面との離間距離を計算する段階と、該離間距離に基づいて前記内周面と前記外周面とが互い接近する向きの引力を与える段階とを含むことを特徴とする。【選択図】  図１５
2005351770.tif 000002
【特許請求の範囲】
【請求項１】
マトリックスゴムとフィラーとを含むゴム材料の変形をシミュレーションするゴム材料のシミュレーション方法であって、
少なくとも一つの前記フィラーを有限個の要素に分割したフィラーモデルと、
前記マトリックスゴムを、初期状態で前記フィラーモデルの外周面と一致ししかも前記外周面から離間可能な内周面を有した有限個の要素に分割したマトリックスモデルとを含むゴム材料モデルを設定するステップ、
前記ゴム材料モデルに引張変形の条件を与えて変形計算を行うステップ、及び  前記ゴム材料モデルから必要な物理量を取得するステップを含み、
かつ、前記変形計算を行うステップは、前記マトリックスモデルの内周面と、前記フィラーモデルの外周面との離間距離を計算する段階と、
該離間距離に基づいて前記内周面と前記外周面とが互い接近する向きの引力を与える段階とを含むことを特徴とするゴム材料のシミュレーション方法。
【請求項２】
前記フィラーは、カーボンブラックである請求項１記載のゴム材料のシミュレーション方法。
【請求項３】
前記引力は、ゴム材料モデルの負荷変形時において、前記離間距離の増大に伴って徐々に増大しかつピークに至る漸増領域と、
さらなる前記離間距離の増大によって滑らかに減少しかつ前記離間距離が予め定められた特徴距離で零に至る漸減領域とを含む放物状曲線を描く関数に基づいて決定されることを特徴とする請求項１又は２に記載のゴム材料のシミュレーション方法。
【請求項４】
前記引力は、前記漸減領域から初期状態への除荷変形において、除荷開始時の値から零まで線形に減少することにより、前記ゴム材料モデルの変形計算においてエネルギーロスが導入されることを特徴とする請求項３に記載のゴム材料のシミュレーション方法。
【請求項５】
前記ゴム材料モデルは、前記フィラーモデルの前記外周面と、前記マトリックスモデルの前記内周面との間に定義された界面モデルを含むとともに、
前記界面モデルは、初期状態で厚さが零、かつ、該初期状態で互いに重なる前記外周面と前記内周面との節点が離間することにより該節点が囲む大きさに変形する要素からなることを特徴とする請求項１乃至４のいずれかに記載のゴム材料のシミュレーション方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、ゴム材料のフィラー間の相互作用を解析するのに役立つゴム材料のシミュレーション方法に関する。
【背景技術】
【０００２】
タイヤ、スポーツ用品、その他各種の工業製品に使用されるゴム材料は、負荷を受けると大きく変形し、負荷を完全に取り除くと元の状態に復元する。またゴム材料は、変形に際して、応力と歪とが比例せず、かつ、荷重の負荷時と除荷時とでは応力−歪曲線の経路を異にするループを描く非線形性を持つ。試作の手間及びコストを減じるために、ゴム材料の変形過程をコンピュータによりシミュレーションすることが行われている。従来、ゴム材料のシミュレーション方法としては、下記の特許文献１や非特許文献１等が知られている。
【０００３】
【特許文献１】特開２００２−３６５２０５号公報
【非特許文献１】Ellen M. Arruda and Marry C. Boyce著「 A THREE-DIMENSIONAL CONSTITUTIVE MODEL FOR THE LARGE STRECH BEHAVIOR OF RUBBER ELASTIC MATERIALSS」 Journal of the Mechani cs and Physics of Solids Volume 41, Issue 2, Pages 389-412 (February 1993)
【０００４】
特許文献１は、ゴム材料（粘弾性材料）が歪速度に応じて異なった変形をなす点に着目する。具体的には、予め粘弾性材料の実使用状態を想定した測定条件で、当該粘弾性材料に生じる歪、歪速度、応力などを測定し、擬似的な縦弾性係数（正しくはゴムはフックの法則を満たす弾性体ではない。）と歪ないし歪速度との対応関係を導出する工程が行われる。そして、解析対象である粘弾性材料の製品モデルに対して、所定の歪速度を与えるとともに、上記対応関係から縦弾性係数を適宜計算して変形シミュレーションが行われる。特許文献１の方法は、粘弾性材料の擬似的な縦弾性係数が、歪速度に応じて変化するという試みを含むため、ヒステリシスロスをシミュレーションの中に取り込み得る点で評価できる。
【０００５】
またArrudaらが提案した非特許文献１は、分子鎖網目理論を用いることにより、ゴム材料を高分子レベルにまでモデル化して変形計算することが記載されている。
【発明の開示】
【発明が解決しようとする課題】
【０００６】
工業製品として使用されているゴム材料には、通常、マトリックスゴムの中にカーボンブラック等のフィラー（充填剤）が配合される。ゴムに対するフィラーの配合量が増やすと、エネルギーロスが大きくなる。図２４には、発明者らによる実験結果の一例が示される。図において、カーボンブラックの配合量ＣＢが増えると、エネルギーロスに相当するヒステリシスループの面積が大きくなることが分かる。
【０００７】
図２５には、フィラーとしてカーボンブラックが配合されたゴム材料に引張変形を与えたときのマトリックスゴムとカーボンブラックとの界面部分の顕微鏡拡大写真が示される（「ゴムのカーボンブラック補強  その３．カーボン補強に関する最近の重要報告」深堀美英  日本ゴム協会誌  第７７巻  第３号（２００４）  Ｐ．１０７より引用）。この写真において黒い球状のものがカーボンブラックである。やや薄い黒色部分は、マトリックスゴムである。引張変形時、マトリックスゴムの一部は、カーボンブラックとの界面から剥離することが知られている。これにより、ゴム材料の中に、白抜き部分で示される空隙（void）が生成される。このような空隙の形成に際して、マトリックスゴムとフィラーとの界面部分に滑りが生じ、そこではエネルギーロスが生じていると考えられている。また除荷に伴って空隙は消失するが、この際にも、分子鎖の結合によってエネルギーロスが生じると考えられている。フィラーを多量に配合されたゴム材料は、マトリックスゴムとフィラーとの界面部分が増える。この考えは、フィラーを多く配合されたゴムはエネルギーロスが大きくなるという実験結果に符号する。
【０００８】
然るに、従来のゴム材料のシミュレーションでは、このようなフィラーとマトリックスゴムとの界面での剥離を変形計算の中に取り込んだものは存在していない。例えば、本件出願人が提案した特願２００３−３５８１６７号では、マトリックスゴムをモデル化したマトリックスモデルと、フィラーをモデル化したフィラーモデルと、マトリックスモデルとフィラーモデルとの間に界面モデルとを設定し、この界面モデルをマトリックスゴムよりも良く伸びるようにパラメータを調節するものが提案されている。このモデルは、界面を考慮に入れている点で評価しうるが、マトリックスゴムとフィラーとの界面部分における剥離を表現しない。
【０００９】
本発明は、以上のような問題点に鑑み案出なされたもので、少なくとも一つの前記フィラー粒子を有限個の要素に分割したフィラーモデルと、マトリックスゴムを、前記フィラーモデルを取り囲むとともに前記フィラーモデルの外周面と一致ししかも前記外周面から離間可能な内周面を有する界面要素を含む有限個の要素に分割したマトリックスモデルとを含むゴム材料モデルを設定し、かつその変形計算において、引張変形時の前記界面要素の内周面と、前記フィラーモデルの外周面との離間距離を計算する段階と、該離間距離に基づいて界面モデルの内周面にフィラーモデルの外周面に近づく向きの引力を与える段階とを含ませることを基本として、フィラーとマトリックスゴムとの界面における剥離を表現し実際のゴムの変形挙動に近いシミュレーションを可能としうるゴム材料のシミュレーション方法を提供することを目的としている。
【課題を解決するための手段】
【００１０】
本発明のうち請求項１記載の発明は、マトリックスゴムとフィラーとを含むゴム材料の変形をシミュレーションするゴム材料のシミュレーション方法であって、少なくとも一つの前記フィラーを有限個の要素に分割したフィラーモデルと、前記マトリックスゴムを、初期状態で前記フィラーモデルの外周面と一致ししかも前記外周面から離間可能な内周面を有した有限個の要素に分割したマトリックスモデルとを含むゴム材料モデルを設定するステップ、前記ゴム材料モデルに引張変形の条件を与えて変形計算を行うステップ、及び前記ゴム材料モデルから必要な物理量を取得するステップを含み、かつ、前記変形計算を行うステップは、前記マトリックスモデルの内周面と、前記フィラーモデルの外周面との離間距離を計算する段階と、該離間距離に基づいて前記内周面と前記外周面とが互い接近する向きの引力を与える段階とを含むことを特徴としている。
【００１１】
また請求項２記載の発明は、前記フィラーは、カーボンブラックである請求項１記載のゴム材料のシミュレーション方法である。
【００１２】
また請求項３記載の発明は、前記引力は、ゴム材料モデルの負荷変形時において、前記離間距離の増大に伴って徐々に増大しかつピークに至る漸増領域と、さらなる前記離間距離の増大によって滑らかに減少しかつ前記離間距離が予め定められた特徴距離で零に至る漸減領域とを含む放物状曲線を描く関数に基づいて決定されることを特徴とする請求項１又は２に記載のゴム材料のシミュレーション方法である。
【００１３】
また請求項４記載の発明は、前記引力は、前記漸減領域から初期状態への除荷変形において、除荷開始時の値から零まで線形に減少することにより、前記ゴム材料モデルの変形計算においてエネルギーロスが導入されることを特徴とする請求項３に記載のゴム材料のシミュレーション方法である。
【００１４】
また請求項５記載の発明は、前記ゴム材料モデルは、前記フィラーモデルの前記外周面と、前記マトリックスモデルの前記内周面との間に定義された界面モデルを含むとともに、前記界面モデルは、初期状態で厚さが零、かつ、該初期状態で互いに重なる前記外周面と前記内周面との節点が離間することにより該節点が囲む大きさに変形する要素からなることを特徴とする請求項１乃至４のいずれかに記載のゴム材料のシミュレーション方法である。
【発明の効果】
【００１５】
本発明のゴム材料のシミュレーション方法では、変形計算において、マトリックスモデルの内周面と、フィラーモデルの外周面との離間距離を計算する段階と、該離間距離に基づいて前記内周面にフィラーモデルの外周面に近づく向きの引力を与える段階とを含む。したがって、引張変形時、マトリックスモデルがフィラーモデルから剥離する現象を変形計算の中に取り込むことができる。これは、実際のゴム変形により近いシミュレーションを可能とするのに役立つ。またマトリックスモデルの内周面には、フィラーモデルの外周面の離間距離に基づいて定められる引力が与えられる。このため、フィラーモデルとマトリックスゴムとを界面で互いに関連づけ、その相互作用についても種々設定することができる。
【００１６】
また請求項４記載の発明のように、引力を、漸減領域から初期状態への除荷変形において、引力を線形に減少させる場合、エネルギーロスを変形計算に導入できる。したがって、より実際のゴム材料の変形挙動に近いシミュレーションが可能になる。
【００１７】
また上述のような変形計算は、請求項５に記載されるように、ゴム材料モデルに、フィラーモデルの外周面と、前記マトリックスモデルの内周面との間に定義されかつ初期状態で厚さが零、かつ、該初期状態で互いに重なる前記外周面と前記内周面との節点が変位することにより該節点が囲む大きさに変形する要素からなる界面モデルを含ませることにより、容易になし得る。
【発明を実施するための最良の形態】
【００１８】
以下、本発明の実施の一形態を図面に基づき説明する。
図１には、本発明のシミュレーション方法を実施するためのコンピュータ装置１が示される。このコンピュータ装置１は、本体１ａと、入力手段としてのキーボード１ｂ、マウス１ｃと、出力手段としてのディスプレイ装置１ｄとを含む。本体１ａには、図示していないが、演算処理装置（ＣＰＵ）、ＲＯＭ、作業用メモリー、磁気ディスクなどの大容量記憶装置、ＣＤ−ＲＯＭやフレキシブルディスクのドライブ１ａ１、１ａ２が設けられている。そして、前記大容量記憶装置には後述する本発明のシミュレーション方法を実行するための処理手順（プログラム）が記憶されている。したがって、このコンピュータ装置１は、ゴム材料のシミュレーション装置を構成する。
【００１９】
図２には、シミュレーション方法の処理手順の一例が示される。本実施形態では、先ずフィラーが配合されたゴム材料モデルが設定される（ステップＳ１）。図３には、微視構造（ユニットセル）としてのゴム材料モデル２の初期状態が視覚化して示され、図４には、さらにそのＡ部が拡大して示されている。
【００２０】
該ゴム材料モデル２は、解析しようとするゴム材料（実在するか否かを問わない）の微小領域が、有限個の小さな要素（例えば２ａ、２ｂ、２ｃ…）を用いて分割されたものである。各要素２ａ、２ｂ、２ｃ…は、数値解析が可能に定義される。数値解析が可能とは、例えば有限要素法、有限体積法、差分法又は境界要素法といった数値解析法により、各要素ないし系全体についての変形計算が可能な決まりに従って必要なパラメータが設定されていることを意味する。具体的には、各要素２ａ、２ｂ、２ｃ…について、節点の座標値、要素の形状、材料特性値などが定義される。各要素２ａ、２ｂ、２ｃ…には、例えば２次元平面としての三角形ないし四辺形の要素、３次元要素としては、例えば４ないし６面体の要素が好ましく用いられる。このように、ゴム材料モデル２は、前記コンピュータ装置１にて取り扱い可能な数値データを構成する。
【００２１】
また本実施形態ではゴム材料モデル２の変形計算を平面歪状態で行うものが例示される。したがってＺ方向（紙面と垂直方向）には、歪が生じない。またこの例ではユニットセルとしてのゴム材料モデル２の一辺の長さは、縦横それぞれ３５nm×３５nmの正方形としている。
【００２２】
この実施形態のゴム材料モデル２は、少なくとも一つのフィラーを有限個の要素に分割してモデル化されたフィラーモデル３と、マトリックスゴムを複数個の要素で分割してモデル化されたマトリックスモデル４とを含む。
【００２３】
図５には、フィラーの一例としてゴム中に充填されたカーボンブラックの電子顕微鏡写真を模写した略図が示される。カーボンブラック６は、炭素原子からなる直径数１０nm程度のほぼ球状ないし卵形状の一次粒子７を最小構成単位とし、これが単体又は図のように不規則に３次元的に結合した葡萄の房状の構造（これは、「二次粒子」である）を持つ。本実施形態のフィラーモデル３は、フィラーとして、カーボンブラックの一次粒子を球状にモデル化したものが例示される。
【００２４】
この実施形態のフィラーモデル３は、四辺形の平面要素に分割されモデル化される。四辺形には、四辺の長さが等しい要素及び四辺の長さが等しくない要素の双方が含まれる。またこの実施形態では、変形計算を簡略化するために、図３のｘ軸及びｙ軸に対して線対称にモデル化されている。またフィラーモデル３は、図４に示されるように、外周面Ｅａを持っている。フィラーモデル３の外周面Ｅａは、該フィラーモデル３の最外側に配置された要素３ａ、３ｂ、３ｃ…の外側に現れる節点３Ａ、３Ｂ、３Ｃ…間をそれぞれ継いでいる直線の辺によって構成される。
【００２５】
また、カーボンブラックに代表される各種のフィラーは、マトリックスゴムに比べるとその数百倍の硬さを持つ。このため、本実施形態ではフィラーモデル３を弾性体として取り扱う。それゆえ、フィラーモデル３は、材料特性値として縦弾性係数が与えられ、変形計算における構成方程式はフックの法則に従うものになる。また、この実施形態では、ユニットセルの中に、フィラーの一次粒子の一つがモデル化されているが、ゴムへのフィラー配合量などを考慮してその個数を増やすこともできる。さらには、フィラーモデル３として、二次粒子の形状をモデル化したものでも良い。またフィラーは、カーボンブラックに限定されるものではなく、例えばシリカ、さらには他のものでも良い。
【００２６】
マトリックスモデル４は、マトリックスゴムを有限個の要素に分割してモデル化されている。この実施形態のマトリックスモデル４は、前記フィラーモデル３を除いた全ての領域に定義される。具体的には、フィラーモデル３の周囲を取り囲み、ゴム材料モデル２の主要部を構成する。この例のマトリックスモデル４は、四辺形の平面要素に分割されてモデル化される。四辺形には、前記同様、四辺の長さが等しい要素、及び四辺の長さが等しくない要素の双方が含まれ得る。ただし、これ以外の他の要素を用いても良いのは言うまでもない。
【００２７】
マトリックスモデル４は、初期状態において、前記フィラーモデル３の外周面Ｅａと一致する内周面Ｅｂを有する。初期状態は、ゴム材料モデル２の体積力を無視し、かつ外力を零とした状態とする。マトリックスモデル４の内周面Ｅｂは、該マトリックスモデル４の最もフィラーモデル３側に配置された要素４ａ、４ｂ、４ｃ…の内側に現れる節点４Ａ、４Ｂ、４Ｃ…間を継ぐ直線の辺によって構成される。初期状態において、マトリックスモデル４の前記各節点４Ａ、４Ｂ、４Ｃ…は、フィラーモデル３の外周面に現れる節点３Ａ、３Ｂ、３Ｃ…と同じ位置にある。したがって、マトリックスモデル４の内周面Ｅｂは、フィラーモデル３の外周面Ｅａと一致し、隙間は零である。
【００２８】
またマトリックスモデル４は、その内周面Ｅｂが、前記フィラーモデル３の外周面Ｅａから離間可能に定義されている。「離間可能に定義」とは、外力の作用により、マトリックスモデル４の内周面Ｅｂが、前記フィラーモデル３の外周面Ｅａから離れることができる余地を持っていれば足りる。マトリックスモデル４の内周面Ｅｂが、フィラーモデル３の外周面Ｅａから離れることにより、その間に計算上の空隙が形成される。フィラーが配合されたゴム材料のシミュレーションにおいて、このような境界条件の設定は、従来には全く見られない新規な構成の一つと言える。
【００２９】
マトリックスモデル４の構成方程式には、本実施形態ではアフィンモデルの分子鎖網目理論が採用される。分子鎖網目理論は、図６（Ａ）、（Ｂ）に示すように、ゴム材料ａは、微視構造として、無秩序に配向された分子鎖ｃが接合点ｂで連結された網目構造を持つとの考えを前提とする。接合点ｂは、例えば分子間の化学的結合であってそれには架橋点などが含まれる。
【００３０】
図６（Ｃ）に示されるように、１本の分子鎖ｃは、複数のセグメントｄから構成される。一つのセグメントｄは、分子鎖網目理論においては繰り返しの最小構成単位である。また一つのセグメントｄは、化学的には図６（Ｄ）に示されるように、炭素原子が共有結合によって連結した複数個のモノマーｆが連結したものと等価である。個々の炭素原子は、原子同士の結合軸の周りで互いに自由に回転しうるため、セグメントｄは全体として曲がりくねるなど様々な形態をとり得る。
【００３１】
分子鎖網目理論では、接合点ｂが原子の揺らぎ周期に対して長時間的には平均位置が変化しないものとし、接合点ｂの回りの摂動を無視する。さらに二つの接合点ｂ、ｂを両端に持つ分子鎖ｃの端−端ベクトル（end-to-end vector ）は、それが埋め込まれているゴム材料の連続体と共に変形するものと仮定される。
【００３２】
上述の非特許文献１において、Aruudaらは、分子鎖網目理論に基づいた８鎖モデルを提案する。この提案は、図７の左側図に示されるように、巨視的にみて、粘弾性材料は微小な８鎖モデルｇが集合した立方体状の網目構造体ｈをなすとの考えを前提とする。一つの８鎖モデルｇは、図７の右側図に拡大して示されるように、分子鎖ｃが立方体の中心に定められた一つの接合点ｂ１から、立方体の各頂点に位置する８つの各接合点ｂ２にそれぞれのびていると仮定される。本実施形態では、マトリックスモデルの構成方程式として、このような８鎖モデルを基本的に採用しこれに修正を加えている。先ずマトリックスモデルの主歪ｂｉと、主ストレッチλｉとの関係は、下記式（１）を満たす。
【００３３】
【数１】
2005351770.tif 000003
【００３４】
ここで、”Ｎ”は前記分子鎖ｃの１本当たりの平均セグメント数、”ＣR ”はマトリックスゴムのせん断モジュラスで、その中で”ｎ”は単位体積当たりの分子鎖ｃの数、”ｋB ”はボルツマン定数、”Ｔ”は絶対温度、”Ｌ”はランジュバン関数である。また速度型の構成方程式は、キルヒホッフの応力速度と歪速度とを用いて下記式（２）で表すことができる。
【００３５】
【数２】
2005351770.tif 000004
【００３６】
また、”σjj”はコーシー応力テンソル、”δij”はクロネッカーのデルタ、Ａijは左ストレッチテンソル、”Κ”はマトリックスモデルの体積を一定に拘束するために導入されるペナルティ法のペナルティパラメータで、本実施形態ではその値を１００とする。
【００３７】
さらに、本実施形態のマトリックスモデル４は、ゴム自身のエネルギーロスを再現するために、Arrudaらの８鎖モデルが修正して用いられる。先ずゴム材料は、複雑に絡み合った前記分子鎖ｃが伸びることによって数百％にも達し得る大きな歪が許容される。言い換えると、荷重負荷における変形過程では、分子鎖ｃの互いに絡み合った部分がほどけて前記接合点ｂの数が減少する。
【００３８】
図８（Ａ）に示されるように、一つの接合点ｂで接合されている分子鎖ｃ１ないしｃ４に矢印方向の引張応力が作用すると、各分子鎖ｃ１ないしｃ４は伸び、接合点ｂは大きな歪を受けて破損（消滅）するものと考える。図８（Ｂ）に示されるように、これまで２本であった分子鎖ｃ１及びｃ２は、あたかも１本の長い分子鎖ｃ５のように振る舞う。分子鎖ｃ３及びｃ４についても同様である。そして、このような現象は、ゴム材料の負荷変形が進むにつれて逐次発生していくものと考えられる。また、その際には多くのエネルギーの損失が生じる。本実施形態では、この様子を前記式（１）に取り入れる。そのために、分子鎖ｃの１本当たりの平均セグメント数Ｎは、負荷変形時と除荷変形時とでは異なった値を示す可変のパラメータとして定義される。ここで、負荷変形時とは、微小時間の間でマトリックスモデル４の歪が増大する変形であり、逆に除荷変形時とは、歪が減少する変形とする。
【００３９】
上記の点について、改めて、図７の巨視的な３次元の網目構造体ｈを参照する。網目構造体ｈは、幅方向、高さ方向及び奥行き方向にそれぞれ８鎖モデルｇがｋ個結合したものとする。ただし、ｋは十分に大きい数とする。該網目構造体ｈに含まれる接合点ｂの総数を「からみ数」として符号”ｍ”で表すと、からみ数ｍは、式（３）で、また網目構造体ｈに含まれる分子鎖ｃの数、即ち、マトリックスモデル４の単位体積中に含まれる分子鎖の数”ｎ”は、式（４）でそれぞれ表し得る。
ｍ＝（ｋ＋１）3   ＋ｋ3     …式（３）
ｎ＝８ｋ3                   …式（４）
【００４０】
ここで、ｋは十分に大きい数であるから、ｋの３次の項以外を省略すると、上式はそれぞれな式（５）及び（６）で表し得る。
ｍ＝２ｋ3                   …式（５）
ｎ＝８ｋ3                   …式（６）
【００４１】
さらに式（５）及び（６）の関係から、からみ数ｍは、ｎを用いて、式（７）で表し得る。
ｍ＝ｎ／４                  …式（７）
【００４２】
さらに、マトリックスモデル４は、変形してもその中に含まれる分子鎖のセグメントの総数ＮA   は変化しないため、式（８）及び（９）が成り立つ。
ＮA ＝ｎ・Ｎ                …式（８）
Ｎ＝ＮA ／ｎ＝ＮA ／４ｍ    …式（９）
【００４３】
式（９）から明らかなように、荷重の負荷過程及びこれに続く歪の回復過程において、前記平均セグメント数Ｎを可変のパラメータとするためには、ゴムの分子鎖ｃのからみ数ｍを変化させれば良い。このように、前記平均セグメント数Ｎを、負荷変形時と除荷変形時とで値の異なる可変のパラメータとすることにより、変形計算の中に、マトリックスモデル４自体のエネルギーロスを取り入れることができる。
【００４４】
前記平均セグメント数Ｎは、様々な方法で決定することができる。一例として、負荷変形時、その歪に関連したパラメータに基づいて増大させることは望ましい。前記歪に関連したパラメータとしては、特に限定されるものではないが、例えば、歪、歪速度又は歪の１次の不変量Ｉ1 などがある。本実施形態では、前記平均セグメント数Ｎを、下記式（１０）で定める。これは、前記平均セグメント数Ｎが、当該マトリックスモデル４の各要素個々において、それぞれ歪の１次の不変量Ｉ1 （より詳しくはその平方根であるパラメータλc ）の関数であることを示す。
【００４５】
【数３】
2005351770.tif 000005
【００４６】
式（１０）は、種々の実験によって定めた一例であり、上記ＡないしＥは、いずれも定数である。これは、ゴム試験片の単純な１軸引張試験などの実測結果から容易に定めることができる。例えば、先ず解析対象となるゴム材料の応力−歪曲線を得る。そして、その荷重除荷時の曲線に沿うように前記ｎ、Ｎを定める。これにより、分子鎖のセグメントの総数ＮA （＝ｎ・Ｎ）が決まる。次に、負荷時、除荷時とも分子鎖ｃの総セグメント数ＮA は同一であるため、荷重負荷時の曲線に整合するよう、各歪における平均セグメント数Ｎを求める（このＮは変化させる。）。そして、決定された負荷時の平均セグメントＮに一致するよう、式（１０）のパラメータＡないしＥを決定する。本実施形態では、Ｎ＝１６を使用し、かつ、負荷終了時のＮが除荷時のそれと等しくなるように設定している。本実施形態では、上記定数を次のように設定している。
Ａ＝＋２．９４９３
Ｂ＝−５．８０２９
Ｃ＝＋５．５２２０
Ｄ＝−１．３５８２
Ｅ＝＋０．１３２５
【００４７】
図９には、前記平均セグメント数ＮとパラメータλC   との関係が示される。歪に関するパラメータであるλC が大きくなると、平均セグメント数Ｎは滑らかに増大する。この例では、パラメータλC の上限は２．５である。後述する変形計算においては、マトリックスモデル４の各要素について、負荷変形時ではパラメータλC が常に計算される。計算されたλC は、式（１０）に代入され、当該要素の当該歪状態における平均セグメント数Ｎが計算される。平均セグメント数Ｎは、適宜のタイミングで計算されて、マトリックスモデル４の構成方程式の中に取り込まれる。なお本実施形態では、マトリックスモデル４の除荷変形時においては平均セグメント数Ｎを一定値とする。
【００４８】
また本実施形態のゴム材料モデル２は、さらに界面モデル５を含むものが例示される。界面モデル５は、フィラーモデル３の外周面Ｅａと、マトリックスモデル４の内周面Ｅｂとの間に定義される。界面モデル３は、変形計算では、図２５に示された空隙を表す。
【００４９】
図１０（Ａ）には、ゴム材料モデル２が部分的に拡大して示されている。フィラーモデル３の外周面Ｅａの各節点３Ｂ、３Ｃ、３Ｄ…と、マトリックスモデル４の内周面Ｅｂとに現れる各節点４Ｂ、４Ｃ、４Ｄ…とが重なっている初期状態では、界面モデル５の厚さは零である（図では見えない）。他方、図１０（Ｂ）に示されるように、ゴム材料モデル２に引張変形が与えられ、前記節点３Ｂ、３Ｃ、３Ｄ…と、節点４Ｂ、４Ｃ、４Ｄ…とが離間すると、界面モデル５は、前記節点により囲まれる大きさに変化する。具体的には、初期状態で重なっていた２つの節点３Ｂ、４Ｂ間を結ぶ辺、３Ｃ、４Ｃ間を結ぶ辺、前記フィラーモデル３の外周面Ｅａを形成する辺及びマトリックスモデル４の内周面Ｅｂを構成する辺で囲まれる四辺形要素５ｂ…などに変形する。
【００５０】
また界面モデル５は、フィラーモデル３の外周面Ｅａと、マトリックスモデル４の内周面Ｅｂとの離間距離に基づき、前記マトリックスモデル４の内周面Ｅｂとフィラーモデル３の外周面とが互いに近づく向きの引力（Traction）を前記各面Ｅａ、Ｅｂ（又はその節点）に与える。界面モデル５は、離間可能に定義されたフィラーモデル３の外周面Ｅａと、マトリックスモデル４の内周面Ｅｂとの間の粘着力（粘着摩擦力）を表現し、これらの界面の相互作用（粘着性）をシミュレーションしうる。金属マトリックスから第２相の粒子が剥離してその界面に空隙が形成される場合、界面に作用する粘着力は、例えばneedleman らの下記の実験が参考になる。この文献は、金属組成物を元に開発されたものであるが、そのオーダなどを適宜かえることで、ゴムマトリックスとフィラーとの剥離及びそれに伴う空隙の生成について同様に適用できる。
A Continuum Model for Void Nucleation by Inclusion Debonding / A. needleman (Journa l of Applied Mechanics SEP. 1987 Vol.54 P.525-531)
【００５１】
図１１には本実施形態で採用されている引力と離間距離との関係がグラフで示されている。このグラフは、上記文献等を参考として定められた。縦軸には引力ｔABが表される。横軸には、離間距離ｕABを特徴距離δで除した無次元量（ｕ／δ）を示す。ここで、引力ｔないし離間距離ｕの添え字Ａ、Ｂは、図１２に示されるように、界面モデル５の前記外周面Ｅａにおける２つの節点１、２の中点Ｂと、初期状態にて該節点１，２と重なっていた界面モデル５のマトリックスモデル４の内周面Ｅｂ側の節点３、４の中点Ａとの離間距離であることを示す。また前記「特徴距離δ」は、図１２の前記中点Ａ、Ｂ間の離間距離ｕABを増大させた場合、その引力ｔABが零となるときの限界の長さである。
【００５２】
本実施形態では、引力ｔABは、離間距離ｕABの増大に伴って徐々に増大しかつピーク（σmax ）に至る漸増領域と、さらなる前記離間距離ｕABの増大によって滑らかに減少し、かつ、前記離間距離ｕABが特徴距離δとなったときに零に至る漸減領域とを含む放物状曲線を描く。即ち、前記引力ｔABは、離間距離ｕABの関数である。引力ｔABがピーク値σmax に達するとき、はじめてマトリックスゴムには分子鎖の滑りが生じるものと考えられている。このピーク値σmax は、界面の強さを表し、この例では特徴距離δのほぼ１／３で生じさせている。また、離間距離ｕABが特徴距離δになったときに、引力ｔABは消失する。これにより、中点Ａ、Ｂ間には完全な滑りが生じる。
【００５３】
特徴距離δと引力ｔABのピーク値である界面の強さσmax との２つのパラメータを種々変更することにより、フィラーモデル３とマトリックスモデル４との界面の剥離に対する強さを種々変えることができる。これは、ゴム材料を構成するフィラー及び／又はマトリックスゴムに合わせて、界面の剥離に対する振る舞いを種々調節可能であることを意味する。
【００５４】
次に前記引力ｔABの関数の定式化について述べる。この例では二次元形を示すが３次元形を用いることもできる。離間距離ｕABの状態にある前記中点Ａ、Ｂ間のポテンシャルエネルギーφ（ｕAB）は、上記needleman の基本式を採用して式（１１）で定め得る。
【００５５】
【数４】
2005351770.tif 000006
【００５６】
また中点Ａ、Ｂ間に作用する引力は、前記ポテンシャルエネルギーφ（ｕAB）を離間距離ｕABで偏微分することにより得られる。これは式（１２）に示される。
【００５７】
【数５】
2005351770.tif 000007
【００５８】
上で述べたとおり、マトリックスゴムとフィラとの界面では、空隙が大きくなると滑りが生じる。また、空隙の消失に際しても、分子鎖の結合によるエネルギー損が生じる。そこで、本実施形態において、前記引力ｔABを、前記漸減領域から初期状態への除荷変形において、除荷開始時の値から零まで線形に減少させる。例えば図１１の漸減領域上における任意の点Ｐｘから除荷が開始される場合、図の鎖線の矢印のように、引力ｔABは、離間距離ｕABの減少に応じて、当該点Ｐｘと原点とを結ぶ直線上の値をとりながら減少する。このため、引力ｔABは、実線矢印及び鎖線矢印で示されるように、負荷及び除荷の１サイクルでエネルギーロスを表すループを描く。したがって、本実施形態のゴム材料モデル２の変形計算には、マトリックスモデル４とフィラーモデル３との界面の剥離によるエネルギーロスが導入され、実際の変形挙動により近いものになる。また、パラメータをより適切に設定することにより、実際のゴム材料と近似したエネルギーロスの計算を行うことが可能になる。
【００５９】
除荷時、線形に減少する引力ｔABは、式（１２）の［１−（ｕAB／δ）］2 の項の（ｕA B ／δ）に点Ｐｘの値が代入された１次式を用いて計算される。また、引力ｔABがピーク値σmax を超えていない場合、引力ｔABは、負荷時と同じルートを通って初期の状態へと戻るように定義される。この場合、可逆変形となりエネルギーロスは生じない。上で述べたとおり、引力ｔABがピーク値σmax を超えていない場合、界面での滑りが開始されないため、エネルギーロスが生じないことにも整合する。界面モデル５を用いて計算された引力ｔは、マトリックスモデル４の内周面Ｅｂ及びフィラーモデル３の外周面Ｅａの各節点に作用するのは言うまでもない。
【００６０】
次に本実施形態のシミュレーション方法では、ゴム材料モデル２に引張変形の条件を与えて変形計算が行われる（ステップＳ２及びＳ３）。
【００６１】
引張変形の条件としては、ゴム材料モデル２に引張歪を与える歪速度や、最大歪量などを含む。また、本実施形態のシミュレーションは、均質化法（漸近展開均質化法）に基づいて行われる。解析対象領域が任意の微視構造の繰り返しによって構成され、その繰り返しの度合いが非常に密な場合、有限要素法で直接領域を離散化できないことがある。均質化法では、図１３に示されるように、ゴム材料モデル２の微視構造（図３で示されたゴム材料モデル２のユニットセル）を縦横に周期的に持っているゴム材料全体Ｍを表現する巨視的スケールｘI と、前記微視構造を表現する微視的スケールｙI との独立した２変数が用いられる。そして、均質化法では、微視的スケールｙI と巨視的スケールｘI という異なる尺度の場におけるそれぞれ独立した変数を漸近展開することにより、図３に示した微視構造のモデル構造を反映させたある程度の大きさを持つゴム材料全体Ｍの平均的な力学応答を求め得る。
【００６２】
漸近展開均質化法自体は、数値計算法の中で既に確立したものであり（例えば下記文献参照）、かつ、本発明に必須の内容ではないため、ここでは詳細に触れない。  Higa,Y.and Tomita,Y,,Computational Prediction of Mechanical Properties of Nickel-based superalloy with gamma Prime Phase Precipitates,Proceedings of ICM8(Victoria,B.C.,Canada),Advance Materials and Modeling of Mechanical Behavior,(Edited by Ellyin,F,and Proven,J.W.),III(1999),1061-1066,Fleming Printing Ltd.. 
比嘉吉一，冨田佳宏，粒子強化型複合材の均質化法による変形挙動のモデル化とシミュレーション，日本機械学会論文集，A66(2000),1441-1446.
【００６３】
また本実施形態では、図１３に示されたなゴム材料モデルＭを１００nm×１００nmの矩形状とし、かつ、巨視的構造において一様な一軸の引張歪Γ2 （その歪速度を１．０×１０-5／ｓとする。）が生じるように、図３に示された微視構造（ユニットセル）に一定の歪速度を与える。また図１３のｘ2 方向に平均歪速度の引張の負荷条件を加え、該ｘ2 方向の歪Γ2 が０．４０に達すると、同じ歪速度で歪を零まで漸減させる除荷条件を定めている。
【００６４】
図１４には、変形計算（変形シミュレーション）のより具体的な処理手順の一例が示される。変形計算では、先ずデータがコンピュータ装置１に入力される（ステップＳ３１）。入力されるデータには、ゴム材料モデル２を構成する数値データや、前記予め設定された変形ないし境界条件等が含まれる。
【００６５】
次に各要素の剛性マトリックスが作成し（ステップＳ３２）、しかる後、全体構造の剛性マトリックスを組み立てる（ステップＳ３３）。全体構造の剛性マトリックスには、既知節点の変位、節点力が導入され（ステップＳ３４）、剛性方程式の解析が行われる。そして、未知節点変位が決定され（ステップＳ３５）、各要素の歪、応力、主応力といった物理量を計算し、出力する（ステップＳ３６ないし３７）。なおステップＳ３８では、計算を終了させるか否かの判定がなされ、否定的である場合には、ステップＳ３２以降を繰り返す。
【００６６】
上述の変形計算は、例えば有限要素法を用いたエンジニアリング系の解析アプリケーションソフトウエア（例えば米国リバモア・ソフトウェア・テクノロジー社で開発されたＬＳ−ＤＹＮＡなど）を用いて行うことができる。
【００６７】
また、本実施形態の変形計算では、定数等を次のように設定している。
【００６８】
＜マトリックスモデル＞
ＣR ＝０．４（ＭＰａ）
Ｎ＝１６
ｎ＝１．０×１０26
Ｔ＝２９６Ｋ
ｋB ＝１．３８０６６×１０-29 
ＮA ＝ｎ・Ｎ＝１．６×１０27
負荷変形時の平均セグメント数Ｎ
Ｎ＝１６
除荷変形時の平均セグメント数Ｎ（一定）
Ｎ＝１６
【００６９】
なおこの例では、負荷変形時の平均セグメント数Ｎを一定としており、界面のエネルギーロスだけに着目しているが、例えば前記平均セグメント数Ｎを下式（定数は適宜変更しうる。）によって定義することでマトリックスモデルのエネルギーロスをも加味することができる。
Ｎ＝-3.2368+20.6175 λc-21.8168 λc 2  +10.8227 λc 3  -1.9003  λc 4 
またゴム材料モデルは、厚さ方向（図３のＺ軸方向）に変化しないように前記３次元８鎖モデルが用いられている。
【００７０】
＜フィラーモデル＞
縦弾性係数Ｅ：１００ＭＰａ
ポアソン比：０．３
【００７１】
＜界面モデル＞
特徴距離δ＝２５nm
界面強さ（ピーク値）σmax ＝０．１ＭＰａ
法線方向ｎとせん断方向ｓとの強さ比γ＝１．０
なお前記特徴距離δや界面強さσmax などは、実際の材料の実験結果などを参考として定めても良いし、またその結果に整合するように適宜設定することもできる。
【００７２】
また、コンピュータ装置１は、前記変形計算の結果から必要な物理量を取得することができる（ステップＳ４）。図１５（Ａ）〜（Ｄ）には、初期状態からマトリックスモデル４とフィラーモデル４との剥離までの変形挙動が可視化して示されている。ここでの物理量は、各要素の歪と応力とを含む。変形計算は、微小な時間きざみで行われ、適当な時間間隔で各要素の歪と応力とを取得することにより、容易にこのような可視化を行うことができる。特に大きな応力が作用している要素には濃いグレーの色彩を付している。
【００７３】
図１５に示されるように、フィラーモデル３とマトリックスモデル４との剥離は、フィラーモデル３の上、下の極で開始され、それは徐々に赤道へと広がる。この空隙の生成により、マトリックスモデル４の応力集中は解放されることが分かる。さらなる引張変形により、空隙（界面モデル）は引張軸方向に大きく成長し、マトリックスモデル４の内周面Ｅｂに高い伸びと応力とを引き起こす。
【００７４】
また図１６には、他の物理量として、縦軸にゴム材料モデル２の巨視的構造での平均応力Σ2 （下記式（１３）により計算する。）が、横軸にその引張歪Γ2 がそれぞれ示されている。ここでは、引張歪Γ2 の上限は０．４である。なお図中のグラフにおいて、実線は界面が剥離しないように境界条件を設定したものである。それ以外の鎖線のグラフのものは、界面強さσmax を０．１、０．２及び０．３ＭＰａにそれぞれ変化させている。界面でマトリックスモデル４とフィラーモデル３とが剥離しないものでは、応力と歪とがほぼ比例し、エネルギーロスが生じていない。しかし、界面での剥離を取り入れ、しかもそこでエネルギーロスを生じさせた実施例では、いずれもエネルギーロスが生じていることが分かる。
【００７５】
【数６】
2005351770.tif 000008
”Ａ”はユニットセルの面積、ａi は要素ｉの面積、σi 22は要素ｉの引張応力である。
【００７６】
また図１７には、さらに変形計算によって得られた物理量の一例を示し、縦軸にゴム材料モデル２の巨視的構造での対数体積歪Ｅv が、横軸にその引張歪Γ2 がそれぞれ示されている。ここでは、引張歪Γ2 の上限は０．４である。なお図中のグラフにおいて、凡例は、図１６のものと同じである。また太い実線で実際のゴム材料の試験結果が示される。このテストでは、界面強さσmax ＝０．３ＭＰａとしたものが、実際のゴム材料の試験結果に非常に近いことが確認できる。また界面強さを大きくした例では、界面により高い拘束力が働くため、体積歪が非常にゆっくりと増加することがシミュレーションされている。しかし、引力が零となり界面での完全な滑りが生じる段階では、急激な体積歪が生じていることが分かる。
【００７７】
図１８には、界面モデル５の前記特徴距離δを違えた４種類について変形計算を行った結果が示される。左縦軸にはゴム材料モデル２の巨視的構造での平均応力Σ2 、右縦軸にはゴム材料モデルの対数体積歪Ｅv 、及び横軸にはゴム材料モデルの引張歪Γ2 がそれぞれ示されている。ここでは、引張歪Γ2 の上限は０．４である。特徴距離δを大きくした例では、界面での完全な滑りが開始されるまでの距離が大きい。従って、上の界面強さの結果と同様に、体積歪が非常にゆっくりと増加することがシミュレーションされていることが分かる。なお図中ｌ0 は、ゴム材料モデルM の縦の長さ１００nmである。
【００７８】
また図１９には、他の実施結果が示される。
この例では、フィラーモデル３が略楕円形状にモデル化されたものが例示される。また楕円形をなすフィラーモデル３の長径の軸Ｇは、図で示されるｘ軸に対して角度θで傾いている。この例のフィラーモデル３は、例えばカーボンブラックの二次粒子をモデル化したものとして用いられる。
【００７９】
図２０及び図２１には、フィラーモデル３の短径ａと、長径ｂとの比（Ｌａ／Ｌｂ）で表されるアスペクト比を変化させて、上記と同様の条件で変形計算を行った結果が示される。この例では長径の軸がなす角度θは０゜である。図２０では、縦軸に巨視構造の平均応力Σ2 が、横軸にその引張歪Γ2 がそれぞれ示されている。図２１では、縦軸が対数体積歪Ｅv である。いずれも引張歪Γ2 の上限は０．４である。図２０ないし２１の結果から明らかなように、高いアスペクト比のフィラーモデルを含むゴム材料モデルでは、全ての歪区間において、引張変形に対する低い抵抗性と、急速な体積歪の増加が見られる。
【００８０】
さらに、図２２及び図２３には、フィラーモデル３の短径Ｌａと、長径Ｌｂとの比（Ｌａ／Ｌｂ）で表されるアスペクト比を２．０で一定とし、その長径の角度θを変化（θ＝０〜９０゜）させて、上記と同様の条件で変形計算を行った結果が示される。図２２では、縦軸に巨視構造の平均応力Σ2 が、横軸にその引張歪Γ2 がそれぞれ示されている。図２３では、縦軸が対数体積歪Ｅv である。これらの結果から明らかなように、フィラーモデル３の長径の角度θが増すと、全ての歪区間において、引張変形に対する高い抵抗性が見られ、かつ、エネルギーロスも小さくなる。また角度θが大きいと、体積歪のゆっくりとした増加が引き起こされている。
【産業上の利用可能性】
【００８１】
上述のようなゴム材料のシミュレーションは、例えばゴム材へのフィラーの配合割合、フィラーの種類を種々違えたゴム材について、エネルギーロスを含めた粘弾性特性をより詳しく評価することができる。具体的な用途として、例えばフィラーモデルの粒子の大きさやアスペクト比などを種々異ならせてシミュレーションを行い、最適なエネルギーロスが得られる粒子形状を決定し、この形状に基づいて実際のフィラーを試作し、評価することができる。また物性を異ならせたフィラーとマトリックスゴムとを種々組み合わせてシミュレーションを行い、最適なエネルギーロスが得られる配合を調べることもできる。また本発明のシミュレーションでは、フィラーモデルとマトリックスモデルとの剥離というより実際の挙動に近いシミュレーションが可能であり、その結果を利用して実際のゴム材料の試作に先立つ解析を高い精度で行うことも可能とする。
【図面の簡単な説明】
【００８２】
【図１】本実施形態で用いたコンピュータ装置の一例を示す斜視図である。
【図２】本実施形態の処理手順を示すフローチャートである。
【図３】ゴム材料モデル（微視構造）の一実施形態を示す線図である。
【図４】そのＡ部拡大図である。
【図５】カーボンブラックの形状を示す線図である。
【図６】（Ａ）は粘弾性材料、（Ｂ）はその分子鎖１構造を説明する線図、（Ｃ）は１本の分子鎖の拡大図、（Ｄ）はセグメントの拡大図である。
【図７】左側図は粘弾性材料の網目構造体を示す斜視図、右側図は８鎖モデルの一例を示す斜視図である。
【図８】（Ａ）、（Ｂ）は分子鎖の接合点の破断を説明する線図である。
【図９】パラメータλc と分子鎖１本当たりの平均セグメント数Ｎとの関係を示すグラフである。
【図１０】（Ａ）は初期状態、（Ｂ）は変形後のゴム材料モデルの界面の状態を示す略図である。
【図１１】引力と離間距離との関係を示すグラフである。
【図１２】界面の拡大説明図である。
【図１３】均質化法を説明する微視構造と全体構造との関係を示す概念図である。
【図１４】変形計算の手順を示すフローチャートである。
【図１５】（Ａ）〜（Ｄ）は変形計算の結果を視覚的に示すゴム材料モデルの変形図である。
【図１６】ゴム材料モデルの変形計算結果として平均応力と引張歪との関係を示すグラフである。
【図１７】ゴム材料モデルの変形計算結果として体積歪と引張歪との関係を示すグラフである。
【図１８】特徴距離δを違えた場合のゴム材料モデルの変形計算結果として、平均応力、体積歪及び引張歪との関係を示すグラフである。
【図１９】ゴム材料モデルの他の実施形態を示す線図である。
【図２０】フィラーモデルのアスペクト比を異ならせた場合のゴム材料モデルの変形計算結果として平均応力と引張歪との関係を示すグラフである。
【図２１】フィラーモデルのアスペクト比を異ならせた場合のゴム材料モデルの変形計算結果として体積歪と引張歪との関係を示すグラフである。
【図２２】フィラーモデルの長径の角度を異ならせた場合のゴム材料モデルの変形計算結果として平均応力と引張歪との関係を示すグラフである。
【図２３】フィラーモデルの長径の角度を異ならせた場合のゴム材料モデルの変形計算結果として体積歪と引張歪との関係を示すグラフである。
【図２４】カーボンブラックの配合量を違えたゴム材料の引張試験結果を示すグラフである。
【図２５】カーボンブラック配合ゴムの顕微鏡写真である。
【符号の説明】
【００８３】
１  コンピュータ装置
２  ゴム材料モデル
３  マトリックスモデル
４  フィラーモデル
５  界面モデル
【図１】
2005351770.tif 000009
【図２】
2005351770.tif 000010
【図３】
2005351770.tif 000011
【図４】
2005351770.tif 000012
【図５】
2005351770.tif 000013
【図６】
2005351770.tif 000014
【図７】
2005351770.tif 000015
【図８】
2005351770.tif 000016
【図９】
2005351770.tif 000017
【図１０】
2005351770.tif 000018
【図１１】
2005351770.tif 000019
【図１２】
2005351770.tif 000020
【図１３】
2005351770.tif 000021
【図１４】
2005351770.tif 000022
【図１５】
2005351770.tif 000023
【図１６】
2005351770.tif 000024
【図１７】
2005351770.tif 000025
【図１８】
2005351770.tif 000026
【図１９】
2005351770.tif 000027
【図２０】
2005351770.tif 000028
【図２１】
2005351770.tif 000029
【図２２】
2005351770.tif 000030
【図２３】
2005351770.tif 000031
【図２４】
2005351770.tif 000032
【図２５】
2005351770.tif 000033