(57)【要約】【課題】計算時間を飛躍的に短縮できる個別要素法の計算方法を提供すること。【解決手段】流体中の解析対象物を離隔や接触を繰り返す複数の小要素の集合体として取り扱う個別要素法の計算方法において、小要素間の抵抗を示す抗力係数と、流体中の小要素の慣性力を示す慣性力係数と、小要素の外部との摩擦抵抗を示す摩擦係数と、流体の流速ベクトルとを設定変数として抽出し、前記設定変数に確率分布を与え、小要素の変位ベクトルをテーラー展開し、前記設定変数とテーラー展開した前記変位ベクトルを運動方程式に代入して摂動展開をおこなって変形量の期待値及び標準偏差を算出する方法である。【選択図】        図２
2004310381.tif 000002
【特許請求の範囲】
【請求項１】
解析対象物を複数の小要素の非連続な集合体として取り扱う個別要素法の計算方法において、
計算に必要な設定変数を抽出し、前記設定変数の一部又は全部に確率分布を与え、
前記設定変数を運動方程式に代入して摂動展開をおこなって解析値の統計量を算出する、
個別要素法の計算方法。
【請求項２】
解析対象物を複数の小要素の非連続な集合体として取り扱う個別要素法の計算方法において、
計算に必要な設定変数を抽出し、前記設定変数の一部又は全部に確率分布を与え、
小要素の変位ベクトルをテーラー展開し、
前記設定変数とテーラー展開した前記変位ベクトルを運動方程式に代入して摂動展開をおこなって変形量の期待値及び標準偏差を算出する、
個別要素法の計算方法。
【請求項３】
解析対象物を離隔や接触を繰り返す複数の小要素の集合体として取り扱う個別要素法の計算方法において、
小要素間の減衰を示すダッシュポットと、小要素間の反力を示すバネ定数と、小要素の外部との摩擦抵抗を示す摩擦係数と、外力とを設定変数として抽出し、前記設定変数の一部又は全部に確率分布を与え、
小要素の変位ベクトルをテーラー展開し、
前記設定変数とテーラー展開した前記変位ベクトルを運動方程式に代入して摂動展開をおこなって変形量の期待値及び標準偏差を算出する、
個別要素法の計算方法。
【請求項４】
流体中の解析対象物を離隔や接触を繰り返す複数の小要素の集合体として取り扱う個別要素法の計算方法において、
小要素間の抵抗を示す抗力係数と、流体中の小要素の慣性力を示す慣性力係数と、小要素の外部との摩擦抵抗を示す摩擦係数と、流体の流速ベクトルとを設定変数として抽出し、前記設定変数に確率分布を与え、
小要素の変位ベクトルをテーラー展開し、
前記設定変数とテーラー展開した前記変位ベクトルを運動方程式に代入して摂動展開をおこなって変形量の期待値及び標準偏差を算出する、
個別要素法の計算方法。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
本発明は、沿岸・海洋構造物の捨石構造物の挙動や亀裂を伴う岩盤の挙動の解析などに使用される個別要素法の計算方法に関するものである。
【０００２】
【従来の技術】
例えば、護岸の捨石構造物は、捨石を積み上げて構築するだけであって、捨石間を連結させることはない。このため、波浪などの外力が捨石に作用した場合、捨石構造物はコンクリート壁のように連続体として挙動するのではなく、個々の捨石が離隔したり再接触したりしながらそれぞれ別の動きをすることになる。このような複雑な挙動を示す構造物は、従来、模型実験によって変形量や被災の不確定性を検討していたが、多大なコストがかかるという問題があった。
【０００３】
そこで、このような検討を数値解析で行うことが提案され、その解析方法のひとつに個別要素法がある。個別要素法は、解析対象物の媒質を複数の小要素の集合体として取り扱う非連続体解析法である。例えば捨石構造物を例にして説明すると、ひとつの捨石にひとつの要素を対応させ、要素間に配置したバネとダッシュポットで要素間の力の作用を算定し、要素の運動方程式に逐次代入していくことで解析をおこなう（図２参照）。
一般的な数値解析では、解析に必要な設定変数（パラメタ）に対して解析結果が算出され、設定変数の設定値と解析結果は１対１に対応する、いわゆる決定論的手法である。そこで、設定変数に様々な値を設定する計算を、破壊確率が１／１００のオーダーの場合は１，０００回繰り返し、破壊確率が１／１０００のオーダーの場合は１０，０００回繰り返す必要がある。そして、蓄積された結果に対して期待値、標準偏差あるいは破壊確率や被災率を算出する。これがいわゆるモンテカルロシミュレーションであり、構造物の変形量や被災を確率的に検討する方法である。
【０００４】
【非特許文献１】
伊藤一教、東江隆夫、勝井秀博、“ＤＥＭ法とＶＯＦ法を用いた粒状体構造物の破壊シミュレーション”、海岸工学論文集、社団法人土木学会、２０００年１０月、第４７巻、ｐｐ．７４６−７５０
【非特許文献２】
伊藤一教、樋口雄一、東江隆夫、勝井秀博、“個別要素法に基づく捨石のランダム性を考慮した潜堤の変形予測手法”、海岸工学論文集、社団法人土木学会、２００１年１０月、第４８巻、ｐｐ．８０６−８１０
【０００５】
【発明が解決しようとする課題】
前記した従来の個別要素法の計算方法にあっては、次のような問題点がある。＜イ＞モンテカルロシミュレーションが適用できるのは、１，０００〜１０，０００回の繰り返し計算を実施しても解析に要する時間が許容できる場合であり、シンプルな問題に限られる。したがって、１回の計算に長時間の解析時間を要する問題に対しては、モンテカルロシミュレーションは現実的でない。
＜ロ＞実際の構造物又はそれを構成する小要素は、様々な形状をしている。また、自然現象に起因する外力や、要素間に設定するバネやダッシュポットなどの設定変数は、ばらつきを有する不確定なものが多い。このため、実際の構造物の解析は複雑な問題となり、上記した理由により数値解析をおこなうには限界があった。
＜ハ＞解析時間を短縮するために、２次元断面を抽出したり、検討モデルを単純化したりする手法が従来から実施されていたが、解析結果まで概算となり、定性的な結果しか得られない場合が多い。
【０００６】
【発明の目的】
本発明は上記したような従来の問題を解決するためになされたもので、計算時間を飛躍的に短縮できる個別要素法の計算方法を提供することを目的とする。
また、実際の構造及び現象に近い状態での解析が容易に実施できる個別要素法の計算方法を提供することを目的とする。
さらに、精度の高い定量的な解析結果を短時間で得ることが可能な個別要素法の計算方法を提供することを目的とする。
本発明は、これらの目的の少なくとも一つを達成するものである。
【０００７】
【課題を解決するための手段】
上記のような目的を達成するために、本発明の個別要素法の計算方法は、解析対象物を複数の小要素の非連続な集合体として取り扱う個別要素法の計算方法において、計算に必要な設定変数を抽出し、前記設定変数の一部又は全部に確率分布を与え、前記設定変数を運動方程式に代入して摂動展開をおこなって解析値の統計量を算出する方法である。
また、解析対象物を複数の小要素の非連続な集合体として取り扱う個別要素法の計算方法において、計算に必要な設定変数を抽出し、前記設定変数の一部又は全部に確率分布を与え、小要素の変位ベクトルをテーラー展開し、前記設定変数とテーラー展開した前記変位ベクトルを運動方程式に代入して摂動展開をおこなって変形量の期待値及び標準偏差を算出する方法である。
【０００８】
さらに、解析対象物を離隔や接触を繰り返す複数の小要素の集合体として取り扱う個別要素法の計算方法において、小要素間の減衰を示すダッシュポットと、小要素間の反力を示すバネ定数と、小要素の外部との摩擦抵抗を示す摩擦係数と、外力とを設定変数として抽出し、前記設定変数の一部又は全部に確率分布を与え、小要素の変位ベクトルをテーラー展開し、前記設定変数とテーラー展開した前記変位ベクトルを運動方程式に代入して摂動展開をおこなって変形量の期待値及び標準偏差を算出する方法である。
【０００９】
そして、流体中の解析対象物を離隔や接触を繰り返す複数の小要素の集合体として取り扱う個別要素法の計算方法において、小要素間の抵抗を示す抗力係数と、流体中の小要素の慣性力を示す慣性力係数と、小要素の外部との摩擦抵抗を示す摩擦係数と、流体の流速ベクトルとを設定変数として抽出し、前記設定変数に確率分布を与え、小要素の変位ベクトルをテーラー展開し、前記設定変数とテーラー展開した前記変位ベクトルを運動方程式に代入して摂動展開をおこなって変形量の期待値及び標準偏差を算出する方法である。
【００１０】
【発明の実施の形態】
以下、図面を参照しながら本発明の実施の形態について説明する。
【００１１】
＜イ＞個別要素法の概要
図１に示すように直径Ｄ（＝２ｒ）、質量ｍの円形要素が振動流下にさらされ、外力として流速の２乗に比例する流体力を考える断面２次元問題を例にとって概略的に説明する。
個別要素法は複数の要素が接触する場合の各要素ｉ，ｊ間に作用する力をバネｋｎ，ｋｓ、ダッシュポッドｃｎ，ｃｓ、摩擦係数μでモデル化し、接触要素に作用する力の総和を評価する（図２参照）。具体的な運動方程式は式（１）、（２）となり、回転運動は式（３）となる。
【００１２】
【数１】
2004310381.tif 000003
【００１３】
ここに、ｋｎ，ｋｓは要素間のバネ定数で添字ｎ，ｓは法線および接線方向を表す。ｃｎ，ｃｓは要素間の減衰定数を表す。ｊは要素ｉと接する要素で、Δｎｉｊ，Δｓｉｊは要素ｉ，ｊの法線および接線方向の相対変位を表し、その微分は時間変化を表す。摩擦係数には静止摩擦、動摩擦および転がり摩擦係数があるが表記を簡便にするため摩擦係数μとして表す。ｌｘ，ｌｚは水平、鉛直成分を表す。Ｉｉ，φｉは要素ｉの慣性モーメントと回転変位で、Ｄｉは減衰定数である。Ｍｉは要素に作用するモーメントである。ｇは重力加速度、ρは流体の密度、Ａは流れに対する要素の投影面積、Ｃｄは抗力係数、ｕとｗは水平・鉛直流速である。
【００１４】
具体的に式（１）、（２）の左辺第２、３項を計算する際には、式（４）〜（１１）の関係式を用いて算出する。
【００１５】
【数２】
2004310381.tif 000004
【００１６】
ここで、ｅｅｎ−ｉｊｌｔ，  ｅｅｓ−ｉｊｌｔは時刻ｔにおける要素ｉ，ｊ間のバネによる作用力を示し、ｅｄｎ−ｉｊｌｔ，  ｅｄｓ−ｉｊｌｔはダッシュポッドによる作用力を示す。ΔＵｎ−ｉｊ，ΔＵｓ−ｉｊおよびΔφｉはΔｔ間の変位増分である。Δｘｉ，Δｚｉは要素ｉのΔｔ間の変位量である。
【００１７】
＜ロ＞本発明における個別要素法の計算方法の概要
本発明の個別要素法の計算方法は、確率個別要素法（Ｓｔｏｃｈａｓｔｉｃ  Ｄｉｓｔｉｎｃｔ  Ｅｌｅｍｅｎｔ  Ｍｅｔｈｏｄ）と呼ぶことができる。本発明の確率個別要素法は、バネ定数、減衰定数、摩擦係数および抗力係数などの設定変数の一部又は全部を確率変数とし、各変数の確率分布が平均値の周りに分布することを仮定する。そして、要素の変位Ｘ，Ｚを微小パラメタにより摂動展開した後、運動方程式に代入して微小パラメタの２次オーダーまでを考慮したＸ，Ｚの摂動展開式を定式化する。ついで、各オーダーの摂動展開式を個別要素法と同様に数値計算し、各オーダーの摂動解を算出する。この摂動解と確率変数の確率分布から決まる４次までのモーメントを用いて、Ｘ，Ｚの期待値と標準偏差を算出する。つまり、従来の個別要素法に対して確率理論における２次近似法を適用したもので、微小パラメタの１次オーダーで打ち切れば、確率理論における１次近似２次モーメント法を個別要素法に適用したことになる。
【００１８】
摂動法を式で表現すると以下のようになる。解析に必要な設定変数（パラメタ）が確率分布を持つ場合、物体の位置もまた確率分布を持つ。ここで，時刻ｔにおける物体の位置Ｘ，Ｚを次式のように示す。
【００１９】
【数３】
2004310381.tif 000005
【００２０】
右辺第１項のＸ０，Ｚ０は解析に必要なパラメタの平均値に対応する物体位置、第２、３項は確率分布の広がりに伴う物体位置を補正する項である。
本発明は、これらの算出結果とパラメタの確率分布によって一義的に決定される確率分布のモーメントを用いて変形や被災に関する統計量（期待値、標準偏差および破壊確率）を算出する方法である。つまり、本発明の特徴は、個別要素法に摂動法と確率理論における１次近似や２次近似法を適用することで、モンテカルロシミュレーションと同等の統計量（期待値、標準偏差および破壊確率）を極めて短時間に算出できる点にある。
【００２１】
＜ハ＞確率個別要素法の定式化
ここでは、ｃやμを固定し抗力係数Ｃｄのみを確率変数として個別要素法の確率的な取り扱いを示す。
まず、Ｃｄは式（１２）のように期待値が０の微小確率変数εで表されるとする。
【００２２】
【数４】
2004310381.tif 000006
【００２３】
確率変数εが確率的に独立で、ｘ，ｚはＣｄの平均値の周りにテーラー展開できると仮定すると、εを摂動パラメタとして式（１３）および（１４）のように摂動展開できる。Ｘｉ（ｍ），Ｚｉ（ｍ）は要素ｉの各オーダーの変位で、添字（）はオーダーを示し、断りがない限り以後同様とする。
【００２４】
【数５】
2004310381.tif 000007
【００２５】
これらを式（１）〜（１１）に代入すると、各式はＸｉ（ｍ），Ｚｉ（ｍ）のべき級数表示ができ、εの各オーダーに対する摂動展開式が得られる。ここでは、εの３乗以上の項を無視する。また、式（８）、（９）および式（１０）はＸｉ（０．．．２），Ｚｉ（０．．．２）で表示されるが、表示を簡略にするため摂動展開式を式（１５）、（１６）のように示す。
【００２６】
【数６】
2004310381.tif 000008
【００２７】
これにより、式（１）の展開式は式（１７）となり、式（２）についても同様に表示できる。式中の（ｍ）はオーダーを示し、Ｆｘ（ｍ）は各オーダーの外力を表す。
【００２８】
【数７】
2004310381.tif 000009
【００２９】
その他の項についてもＦｘ（ｍ）と同様にＸｉ（ｍ），Ｚｉ（ｍ）で表示できる。
式（３）の回転変位については、要素ｉに作用する各オーダーのｅｅｓ−ｉｊｌｔを用いてモーメントを算出し、オーダー毎に回転変位を求める。
【００３０】
さて、確率変数Ｃｄの確率分布に対する１〜４次モーメントが既知であるならば、確率理論における２次近似が適用でき、ｘｉ，ｚｉの期待値Ｅと分散Ｖａｒを算出することができる。式（２２）、（２３）はｘｉに対する期待値と分散の近似式を示す。
【００３１】
【数８】
2004310381.tif 000010
【００３２】
上式から分散を求めれば、標準偏差も同時に求められる。
なお、上記の展開では簡単にするため抗力係数Ｃｄのみを確率変数として説明したが、その他に慣性力係数ＣＭ、摩擦係数μ及び流速ベクトルＵ（ｕ、ｗ）等を確率変数としても、同様に本発明を適用できる。
【００３３】
＜ニ＞確率個別要素法の数値実験による検証
確率個別要素法の近似精度を確かめるため、同一振動流の条件で、Ｃｄの平均値を０．５にして計算した摂動解と、従来の個別要素法によりＣｄを０．２〜１．０の範囲で与えて直接計算した結果（確定解）を比較した。図３に結果を示す。図の横軸εは−０．６の時にＣｄ＝０．２、ε＝１の時にＣｄ＝１．０である。縦軸は要素直径Ｄで無次元化した水平変位ｘ／Ｄである。比較した時刻は移動開始後水平変位が最大になる直前の時刻である。
図よりε＝−０．６や１でも摂動解と確定解の差異は微小である。この結果は、εに微小確率変数を仮定したが、Ｃｄの平均値を０．５にして展開された摂動解がＣｄ＝０．２〜１．０の範囲まで十分な近似精度を有することを示している。また、１次近似と２次近似を比較すると、２次近似の方がわずかに確定解に近いが、その差は有意ではない。これは、Ｃｄが流体力に線形で効くため１次近似でも十分な精度が得られるためである。ここでは２次近似を用いて検証を行った。
つぎに、確率個別要素法で算出される期待値および標準偏差の近似精度を検証するため、図１に示した要素の移動現象を対象に、試行回数１，０００回のモンテカルロシミュレーション結果と比較した。
検証計算に用いたδはＤ／２とＤ／８の２通りとし、Ｕは５ｍ／ｓとした。Ｃｄの確率分布は、図４に示す一様分布、正規確率分布および対数正規確率分布とした。図中には各分布の呼称と統計値を示しており、対数正規確率分布はＬＮ（最頻値、平均値、標準偏差）で表した。ただし、確率個別要素法の計算で用いるＣｄの平均値は、確率分布に関係なく０．６とした。
【００３４】
確率個別要素法（ＳＤＥＭと記す）と従来の個別要素法によるモンテカルロシミュレーション（ＤＥＭｍと記す）で算出した水平変位の経時変化を図５に示す。計算条件はδ＝Ｄ／８、確率分布はＮ３：Ｎ（０．６，０．０８５）とＬ２：ＬＮ（０．５６，０．６３，０．１９）である。図には期待値と期待値±標準偏差を示した。ＳＤＥＭとＤＥＭｍの差異は確率分布Ｎ３、Ｌ５ともに期待値だけでなく標準偏差で示された変動幅も比較的精度良く近似している。
図６は底面条件がδ＝Ｄ／２で、ｔ／Ｔ＝０．６８における期待値、標準偏差および変動係数を比較した図である。横軸は図４に示した確率分布の呼称を示している。期待値は分布形状に関係なくＤＥＭｍがＳＤＥＭを上回っている。これは、ＤＥＭｍの場合要素の移動に伴って変位の分布が変化するが、ＳＤＥＭは変位の分布はほぼ一定であることが起因している。ＳＤＥＭとＤＥＭｍとの差は凹凸間隔δの１〜２倍程度で、この差は要素が移動方向を逆転させる際に凸部の右側で止まるか左側で止まるかによっている。この程度の差は、底面凹凸の分解能に相当し移動量の絶対値が８〜１０δであることからすると誤差は小さいと考えられる。
【００３５】
標準偏差は期待値とは異なり、確率分布による影響が見られる。正規分布Ｎ１〜Ｎ４ではＳＤＥＭとＤＥＭｍとの差異は小さい。Ｌ１〜Ｌ４の対数正規分布では、ＳＤＥＭの計算で用いたＣｄの平均値の０．６と最頻値が近いＬ５やＬ６は差異が小さいが、ＳＤＥＭで用いたＣｄと最頻値が離れるＬ１やＬ４では差異が大きい。また。対数正規分布の場合にはＳＤＥＭのＣｄと最頻値との差だけでなく、分布形が広帯域になることも差異に影響を及ぼす。さらに、一様分布ｕｎｉの場合には、ＳＤＥＭの狭帯近似が成立しないため差異が大きい。
変動係数については標準偏差の結果と類似するが、一様分布ｕｎｉの場合には、期待値と標準偏差の差異が相殺し差異が少なくなった。
対数正規分布で広帯域の場合、平均値と最頻値は差異が大きくなる。ゆえに、ＳＤＥＭの計算で用いるＣｄの平均値に最頻値を用いるほうが適切な場合もあると考えられる。そこで、底面の影響が少ないδ＝Ｄ／８で対数正規分布Ｌ４：ＬＮ（０．７３，０．８５，０．３）の条件に対して、ＳＤＥＭのＣｄの平均値に最頻値０．７３を用いた。その結果、ｔ／Ｔ＝０．６８におけるＤＥＭｍの期待値５．７７ｍに対して５．９１ｍとなり、ほぼ一致するようになった。一方、標準偏差はＤＥＭｍが０．８５６ｍに対して１．１８ｍとなり、若干過大評価となった。したがって、対数正規分布のように平均値と最頻値が異なる場合には、単純に平均値を用いるのではなく、ＳＤＥＭに最頻値を用いることにより、適切な計算結果が得られることがわかる。
【００３６】
ここでＳＤＥＭとＤＥＭｍの計算時間を比較すると同一計算機を用いた１ケース当りの計算時間はＳＤＥＭで１．４分、ＤＥＭｍで８時間であった。つまり、モンテカルロシミュレーションの試行回数Ｎに対して計算量は３／Ｎにまで短縮できる。この特徴は、ＳＤＥＭが変形問題に対する確率的性能照査手法として有力な手法となりうることを示している。
【００３７】
【実施例】
＜イ＞模型実験の結果との比較
図７は混成防波堤（図７のグラフ中の左下）のＡ〜Ｅ区間毎に算出した被災率と波高の関係を示した図で、実験値と本発明の計算値を比較している。このときの周期は１．３ｓである。図７の被災率Ｄｒは次の式で算定できる。
【００３８】
【数９】
2004310381.tif 000011
【００３９】
ここで、ｎは各区間内で被覆石１個分以上移動した個数、Ｎは被覆石の総数である。図７の被覆石総数Ｎは、区間毎の被覆石総数を用いた。
図７（ａ）の区間Ｃにおける被災率をみると、Ｃｄの平均値を１．０とした計算結果が実験結果とよく一致している。ただし、波高１７ｃｍのケースでは破壊確率が過小評価となるが、全体的には比較的よく一致している。また、図７（ｂ）の区間ＢでもＣｄの平均値を１．０とした計算結果は実験結果とよく一致している。ちなみに、実験における変位の許容値Ｘｃは１Ｄである。
ここで実施した計算は、確率変数の平均値、標準偏差および分布形状を便宜的に与えた計算であるため、これらの結果は普遍的なものではない。しかし、ここで用いた確率変数の平均値、標準偏差および確率分布形状で概ね実験結果を再現していることが確認でき、本手法によって確率的な計算結果を再現できることが確認できた。
【００４０】
＜ロ＞従来の計算方法との計算時間の比較
本発明の目的は、モンテカルロ法より短い計算時間で、実験結果に相当する計算結果を得ることにある。そこで、モンテカルロシミュレーションと本発明の確率個別要素法の計算時間についてオーダー比較すると次式のようになる。
【００４１】
【数１０】
2004310381.tif 000012
【００４２】
ここで、ＴＭＯＮ．はモンテカルロシミュレーションによる計算時間で、ＭＭＯＮ．は試行回数である。ＴＳＤＥＭ．は本発明による計算時間で、ＭＳＤＥＭ．は確率変数の数である。式中のＯ（）はオーダーを表す。
一般にＭＭＯＮ．は破壊確率の逆数の１０倍以上を必要とするため、破壊確率が数％の場合にはＭＭＯＮ．は１０００回程度となる。一方、本発明において確率変数を数個設定するならば、計算時間は千分の１〜数百分の１にまで短縮できる。よって、本発明の確率個別要素法は短時間に確率・統計的な結果を解析できる性能照査手法となりうると考えられる。
以上においては、確率変数の分布形状に正規分布を仮定して検討してきた。しかし、確率変数の分布形状やその標準偏差の影響が破壊確率に敏感に現れることがわかっている。ゆえに、確率変数の分布形状やその標準偏差としてどのような形状や値が適切かについては個々の設定変数に対して検討すればよい。
【００４３】
【発明の効果】
本発明の個別要素法の計算方法は以上説明したようになるから次のような効果を得ることができる。
＜イ＞個別要素法に摂動法と確率理論における２次近似法を適用することで、計算時間を飛躍的に短縮できる。
＜ロ＞コストの高い模型実験を行わなくとも、実際の構造及び現象に近い状態での解析が容易に実施できる。
【図面の簡単な説明】
【図１】個別要素法のモデル化の実施例の説明図。
【図２】捨石構造物に個別要素法を適用した場合のモデル化の実施例の説明図。
【図３】本発明の個別要素法の計算方法による計算結果と従来法の計算結果の比較した図。
【図４】検証に用いた確率分布とその統計値を示した図。
【図５】本発明と従来法の水位変位の経時変化を示した図。
【図６】本発明と従来法の期待値及び標準偏差を比較した図。
【図７】被災率に関する実験値と計算値を比較した図。
【図１】
2004310381.tif 000013
【図２】
2004310381.tif 000014
【図３】
2004310381.tif 000015
【図４】
2004310381.tif 000016
【図５】
2004310381.tif 000017
【図６】
2004310381.tif 000018
【図７】
2004310381.tif 000019