(57)【要約】【課題】きしみ割れの発生を定量的に推定することが可能なきしみ割れ推定装置を提供する。【解決手段】車両モデルデータ４０及びレールモデルデータ４１を用いて模擬レール５１上を模擬車両５０が走行した際の模擬レール５１に作用する力を算出する車両運動解析部３１と、車両運動解析部３１により算出された力を用いて模擬レール５１の摩耗の進展速度を算出する摩耗進展解析部３２と、車両運動解析部３１により算出された力を用いて模擬レール５１に生じうるき裂の進展速度を算出するき裂進展解析部３３と、摩耗の進展速度とき裂の進展速度とを用いて模擬レールにきしみ割れが生じるか否かを推定するきしみ割れ推定部３４とを有する。【選択図】図１
2019178944.tif 000002
【特許請求の範囲】
【請求項１】
ゲージコーナー部ならびにレール頭頂部に生じうるきしみ割れの発生を推定するきしみ割れ推定装置であって、
模擬車両のモデルを設定する車両モデルデータ及び模擬レールのモデルを設定するレールモデルデータが格納される記憶部と、前記車両モデルデータ及び前記レールモデルデータを用いて前記模擬レール上を前記模擬車両が走行した際の前記模擬レールに作用する力を算出する車両運動解析部と、前記車両運動解析部により算出された前記力を用いて前記模擬レールの摩耗の進展速度を算出する摩耗進展解析部と、前記車両運動解析部により算出された前記力を用いて前記模擬レールに生じうるき裂の進展速度を算出するき裂進展解析部と、前記摩耗の進展速度と前記き裂の進展速度とを用いて前記模擬レールにきしみ割れが生じるか否かを推定するきしみ割れ推定部とを有することを特徴とするきしみ割れ推定装置。
【請求項２】
前記摩耗進展解析部は、前記車両運動解析部により算出された前記力を用いて前記模擬車両と前記模擬レールとの接触面内における接触応力を算出するレール接触解析部と、前記レール接触解析部により算出された前記接触応力に基づいて前記模擬レールの前記摩耗の進展速度を算出する摩耗量算出部とを有することを特徴とする請求項１に記載のきしみ割れ推定装置。
【請求項３】
前記摩耗進展解析部は、前記摩耗量算出部により算出された摩耗量に基づいて前記レールモデルデータを更新するレールモデルデータ更新部を有し、
前記車両運動解析部は、前記レールモデルデータ更新部により更新された前記レールモデルデータを用いて前記模擬レールに作用する前記力を算出する
ことを特徴とする請求項２に記載のきしみ割れ推定装置。
【請求項４】
前記き裂進展解析部は、前記レール接触解析部により算出された前記接触応力に基づいて前記模擬レールに前記き裂が発生するか否かを判定するき裂発生条件判定部と、前記き裂発生条件判定部により前記き裂が発生すると判定されたら前記模擬レールに生じうる前記き裂の進展速度を算出するき裂進展量算出部とを有することを特徴とする請求項２または３に記載のきしみ割れ推定装置。
【請求項５】
前記き裂進展解析部は、前記き裂発生条件判定部による判定動作に必要なパラメータの初期値を設定する初期値設定部を有し、前記き裂発生条件判定部は前記初期値設定部により設定された前記パラメータの初期値を用いて前記模擬レールに前記き裂が発生するか否かを判定することを特徴とする請求項４に記載のきしみ割れ推定装置。
【請求項６】
模擬車両のモデルを設定する車両モデルデータ及び模擬レールのモデルを設定するレールモデルデータが格納される記憶部を有するきしみ割れ推定装置により実行され、ゲージコーナー部ならびにレール頭頂部に生じうるきしみ割れの発生を推定するきしみ割れ推定方法であって、
前記車両モデルデータ及び前記レールモデルデータを用いて前記模擬レール上を前記模擬車両が走行した際の前記模擬レールに作用する力を算出し、算出した前記力を用いて前記模擬レールの摩耗の進展速度を算出し、算出した前記力を用いて前記模擬レールに生じうるき裂の進展速度を算出し、前記摩耗の進展速度と前記き裂の進展速度とを用いて前記模擬レールにきしみ割れが生じるか否かを推定することを特徴とするきしみ割れ推定方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、レール頭頂部ならびにゲージコーナー部に生じうるきしみ割れの発生を推定するきしみ割れ推定装置及びきしみ割れ推定方法に関するものである。
【背景技術】
【０００２】
レールの車輪との接触部である頭頂部は、車輪からの荷重、摩擦等の影響を受けて金属疲労を受ける。そして、この金属疲労が進行すると、きしみ割れと呼ばれるき裂を生じることが知られている。特に、主に曲線区間の外軌レールのゲージコーナー部において、列車の進行方向に対してある角度をもって狭間隔に並んだきしみ割れと呼ばれるき裂が生じることが知られている。
【０００３】
ゲージコーナー部に発生したきしみ割れは発生位置が等速部に偏っているため、レール探傷車によるき裂の発見も容易でない。また頭側部に発生したきしみ割れは、成長するとき裂同士が連結して剥離に至ることがあるため、レール探傷車による超音波探傷を阻害する恐れがある。従って、レールの保守管理計画を立てる際に、きしみ割れの発生や進展を精度良く予測できれば、適時適切なレール交換等が可能となり、安全性の確保と同時に経済的な効果も期待できる。
【０００４】
本発明者らは、既に、マルチボディダイナミクスによるレール摩耗形状の予測技術及びその解析結果について提案を行った（非特許文献１参照）。また、マルチボディダイナミクス解析ソフトを用いてレールに発生する応力、歪みを求め、求めた応力からレールの頭頂部に生じうるき裂面を探索し、きしみ割れの発生寿命を評価する手法も提案されている（非特許文献２参照）。
【先行技術文献】
【非特許文献】
【０００５】
【非特許文献１】辻江正裕、外２名“ＭＢＤソフトを援用した摩耗形状予測とその考察”、日本機械学会論文集Ｃ編、２０１３年、第７９巻、第８０６号、p.3376−3388
【非特許文献２】松田博之、外５名“レールきしみ割れ発生寿命予測の試み”、第１５回鉄道技術・政策連合シンポジウム予稿集、２００８年、p.537−540
【発明の概要】
【発明が解決しようとする課題】
【０００６】
しかしながら、上述した従来の提案は、それぞれマルチボディダイナミクス解析ソフトを用いてレール摩耗形状の予測を行っているか、あるいは、同様にマルチボディダイナミクス解析ソフトを用いてきしみ割れの発生寿命を評価しており、いずれの提案も、レールと車輪との接触による摩耗・疲労バランスに基づいてきしみ割れの発生を定量的に推定するには至っていない。
【０００７】
そこで、本発明は、きしみ割れの発生を定量的に推定することが可能なきしみ割れ推定装置及びきしみ割れ推定方法を提供することを目的としている。
【課題を解決するための手段】
【０００８】
本発明者は上述の課題に基づき鋭意研究した結果、レール上を車両が走行した際のレールに作用する力を用いてレールの摩耗の進展速度及びレールのき裂の進展速度を算出することができ、そして、これらの関係に基づいてきしみ割れの発生を推定できるとの知見を得るに至った。
【０００９】
そこで本発明のきしみ割れ推定装置は、ゲージコーナー部ならびにレール頭頂部に生じうるきしみ割れの発生を推定するきしみ割れ推定装置であって、模擬車両のモデルを設定する車両モデルデータ及び模擬レールのモデルを設定するレールモデルデータが格納される記憶部と、車両モデルデータ及びレールモデルデータを用いて模擬レール上を模擬車両が走行した際の模擬レールに作用する力を算出する車両運動解析部と、車両運動解析部により算出された力を用いて模擬レールの摩耗の進展速度を算出する摩耗進展解析部と、車両運動解析部により算出された力を用いて模擬レールに生じうるき裂の進展速度を算出するき裂進展解析部と、摩耗の進展速度とき裂の進展速度とを用いて模擬レールにきしみ割れが生じるか否かを推定するきしみ割れ推定部とを有することを特徴とする。
【００１０】
ここで、摩耗進展解析部は、車両運動解析部により算出された力を用いて模擬車両と模擬レールとの接触面内における接触応力を算出するレール接触解析部と、レール接触解析部により算出された接触応力に基づいて模擬レールの摩耗の進展速度を算出する摩耗量算出部とを有する構成とすることができる。
【００１１】
また、摩耗進展解析部は、摩耗量算出部により算出された摩耗量に基づいてレールモデルデータを更新するレールモデルデータ更新部を有し、車両運動解析部は、レールモデルデータ更新部により更新されたレールモデルデータを用いて模擬レールに作用する力を算出する構成とすることができる。
【００１２】
さらに、き裂進展解析部は、レール接触解析部により算出された接触応力に基づいて模擬レールにき裂が発生するか否かを判定するき裂発生条件判定部と、き裂発生条件判定部によりき裂が発生すると判定されたら模擬レールに生じうるき裂の進展速度を算出するき裂進展量算出部とを有する構成とすることができる。
【００１３】
そして、き裂進展解析部は、き裂発生条件判定部による判定動作に必要なパラメータの初期値を設定する初期値設定部を有し、き裂発生条件判定部は初期値設定部により設定されたパラメータの初期値を用いて模擬レールにき裂が発生するか否かを判定する構成とすることができる。
【００１４】
また、本発明のきしみ割れ推定方法は、模擬車両のモデルを設定する車両モデルデータ及び模擬レールのモデルを設定するレールモデルデータが格納される記憶部を有するきしみ割れ推定装置により実行され、ゲージコーナー部ならびにレール頭頂部に生じうるきしみ割れの発生を推定するきしみ割れ推定方法であって、車両モデルデータ及びレールモデルデータを用いて模擬レール上を模擬車両が走行した際の模擬レールに作用する力を算出し、算出した力を用いて模擬レールの摩耗の進展速度を算出し、算出した力を用いて模擬レールに生じうるき裂の進展速度を算出し、摩耗の進展速度とき裂の進展速度とを用いて模擬レールにきしみ割れが生じるか否かを推定することを特徴とする。
【発明の効果】
【００１５】
このように構成された本発明のきしみ割れ推定装置では、車両モデルデータ及びレールモデルデータを用いて模擬レール上を模擬車両が走行した際の模擬レールに作用する力を算出する車両運動解析部と、車両運動解析部により算出された力を用いて模擬レールの摩耗の進展速度を算出する摩耗進展解析部と、車両運動解析部により算出された力を用いて模擬レールに生じうるき裂の進展速度を算出するき裂進展解析部と、摩耗の進展速度とき裂の進展速度とを用いて模擬レールにきしみ割れが生じるか否かを推定するきしみ割れ推定部とを有する。
【００１６】
このようにすることで、模擬レールと模擬車両の車輪との接触状態に基づいてきしみ割れの発生を定量的に推定することができる。
【００１７】
ここで、摩耗進展解析部が、車両運動解析部により算出された力を用いて模擬車両と模擬レールとの接触面内における接触応力を算出するレール接触解析部と、レール接触解析部により算出された接触応力に基づいて模擬レールの摩耗の進展速度を算出する摩耗量算出部とを有するので、摩耗進展解析部による模擬レールの摩耗の進展速度算出処理をより正確に定量的に行うことができ、結果として、きしみ割れの発生をより定量的に推定することができる。
【００１８】
また、摩耗進展解析部が摩耗量算出部により算出された摩耗量に基づいてレールモデルデータを更新するレールモデルデータ更新部を有し、車両運動解析部がレールモデルデータ更新部により更新されたレールモデルデータを用いて模擬レールに作用する力を算出するので、模擬レールの摩耗量を的確にフィードバックして模擬レールに作用する力を算出することができ、結果として、きしみ割れの発生をより定量的に推定することができる。
【００１９】
さらに、き裂進展解析部が、レール接触解析部により算出された接触応力に基づいて模擬レールにき裂が発生するか否かを判定するき裂発生条件判定部と、き裂発生条件判定部によりき裂が発生すると判定されたら模擬レールに生じうるき裂の進展速度を算出するき裂進展量算出部とを有するので、き裂進展解析部による、模擬レールに生じうるき裂の進展速度算出処理をより正確に定量的に行うことができ、結果として、きしみ割れの発生をより定量的に推定することができる。
【００２０】
そして、き裂進展解析部がき裂発生条件判定部による判定動作に必要なパラメータの初期値を設定する初期値設定部を有し、き裂発生条件判定部が初期値設定部により設定されたパラメータの初期値を用いて模擬レールにき裂が発生するか否かを判定するので、模擬レールと模擬車両の車輪との接触状態に基づいてきしみ割れの発生を定量的に推定することができる。
【００２１】
また、本発明のきしみ割れ推定方法は、車両モデルデータ及びレールモデルデータを用いて模擬レール上を模擬車両が走行した際の模擬レールに作用する力を算出し、算出した力を用いて模擬レールの摩耗の進展速度を算出し、算出した力を用いて模擬レールに生じうるき裂の進展速度を算出し、摩耗の進展速度とき裂の進展速度とを用いて模擬レールにきしみ割れが生じるか否かを推定するので、模擬レールと模擬車両の車輪との接触状態に基づいてきしみ割れの発生を定量的に推定することができる。
【図面の簡単な説明】
【００２２】
【図１】本実施の形態であるきしみ割れ推定システムの概略構成を示すブロック図である。
【図２】本実施の形態であるきしみ割れ推定装置の全体動作の一例を示すフローチャートである。
【図３】本実施の形態であるきしみ割れ推定装置のき裂進展解析処理の一例を示すフローチャートである。
【図４】本実施の形態であるきしみ割れ推定装置に用いられる車両モデルデータの一例を示す図である。
【図５】本実施の形態であるきしみ割れ推定装置に用いられるレールモデルデータの一例を示す図である。
【図６】本実施の形態であるきしみ割れ推定装置に用いられるレールモデルデータの一例を示す図である。
【図７】本実施の形態であるきしみ割れ推定装置に用いられるレールモデルデータの一例を示す図である。
【図８】本実施の形態であるきしみ割れ推定装置のレール接触解析部における解析処理を説明するための図である。
【図９】本実施の形態であるきしみ割れ推定装置のレール接触解析部における解析処理を説明するための図である。
【図１０】本実施の形態であるきしみ割れ推定装置のき裂発生条件判定部における判定処理を説明するための図である。
【図１１】本実施の形態であるきしみ割れ推定装置のき裂発生条件判定部における判定処理を説明するための図である。
【図１２】本実施の形態であるきしみ割れ推定装置のき裂発生条件判定部における判定処理を説明するための図である。
【図１３】本実施の形態であるきしみ割れ推定装置のき裂進展量算出部における判定処理を説明するための図である。
【図１４】本実施の形態であるきしみ割れ推定装置のき裂進展量算出部における判定処理を説明するための図である。
【発明を実施するための形態】
【００２３】
（きしみ割れ推定システムの概略説明）
以下、この発明の実施の形態を図面に基づいて説明する。図１は、本実施の形態であるきしみ割れ推定システムＳの概略構成を示すブロック図である。
【００２４】
本実施の形態のきしみ割れ推定システムＳは、きしみ割れ推定装置１０、入力装置１１及び表示装置１２を有する。
【００２５】
本実施の形態のきしみ割れ推定装置（以下、推定装置と省略する）１０は、例えばパーソナルコンピュータ等であり、制御部２０、記憶部２１、入力インタフェース（Ｉ／Ｆ）２２及び出力インタフェース（Ｉ／Ｆ）２３を有する。
【００２６】
制御部２０はＣＰＵ等の演算素子を有する。記憶部２１内に格納されている図略の制御用プログラムが推定装置１０の起動時に実行され、この制御用プログラムに基づいて、制御部２０は記憶部２１等を含む推定装置１０全体の制御を行うとともに、表示制御部３０、車両運動解析部３１、摩耗進展解析部３２、き裂進展解析部３３及びきしみ割れ推定部３４としての機能を実行する。これら各機能部の動作については後述する。
【００２７】
記憶部２１はハードディスクドライブ等の大容量記憶媒体、及びＲＯＭ、ＲＡＭ等の半導体記憶媒体を有する。この記憶部２１には上述の制御用プログラムが格納されているとともに、制御部２０の制御動作時に必要とされる各種データが一時的に格納される。
【００２８】
また、この記憶部２１には、車両モデルデータ４０、レールモデルデータ４１、レール形状データ４２、応力データ４３、摩耗量データ４４、き裂発生初期値データ４５及びき裂進展量データ４６が格納されている。これらデータのうち、応力データ４３、摩耗量データ４４、及びき裂進展量データ４６については少なくとも一時的に記憶部２１に格納されればよい。これらデータの詳細については後述する。
【００２９】
入力インタフェース２２は、推定装置１０に接続された入力装置１１からの各種入力を受け入れ、これを制御部２０に出力する。本実施例の入力装置１１は例えばキーボードやマウス等であり、後述する表示装置１２の表示画面に対して座標指定入力を行いうるものである。
【００３０】
出力インタフェース２３は、制御部２０、特に表示制御部３０から出力された出力信号を受け入れ、これを表示装置１２に出力する。本実施例の表示装置１２は例えば液晶ディスプレイ装置であり、出力インタフェース２３を介して出力された表示制御信号に基づいて図略の表示面に表示画面を表示する。
【００３１】
次に、制御部２０に構成される各機能部の説明をする。
【００３２】
表示制御部３０は、制御部２０及びこの制御部２０に構成される各機能部による処理の結果、表示装置１２に表示画面を生成するための表示制御信号を生成して、この表示制御信号を出力インタフェース２３を介して表示装置１２に出力する。
【００３３】
車両運動解析部３１は、記憶部２１に格納された車両モデルデータ４０及びレールモデルデータ４１を用いて、レールモデルデータ４１により設定された模擬レール上を、同様に車両モデルデータ４０により設定された模擬車両が走行した際の、この模擬レールに作用する力を算出する。具体的には、車両運動解析部３１による算出処理は、ＳＩＭＰＡＣＫやＶＡＭＰＩＲＥ（いずれも商品名）に代表されるマルチボディダイナミクス解析ソフトにより実施可能である。
【００３４】
このマルチボディダイナミクス解析ソフトは、弾性体を含むマルチボディモデルの動的解析を効率よく行うためのツールである。より詳細には、マルチボディダイナミクス解析ソフトは、任意の車両モデルで、車輪やレールの実測形状を適用し、レール通過時の車両の挙動や車輪／レール間の接触を解析する。
【００３５】
本実施の形態である車両運動解析部３１では、対象区間の曲線半径や線形、レール形状等のレールモデルデータ４１及び車両諸元や列車速度等の車両モデルデータ４０をマルチボディダイナミクス解析ソフトに入力し、走行シミュレーション結果から、レールモデルデータ４１により設定された模擬レール上を、同様に車両モデルデータ４０により設定された模擬車両が走行した際の、この模擬レールに作用する力を算出する。本実施の形態である車両運動解析部３１の処理内容、さらにマルチボディダイナミクス解析ソフトそのものは既知のものであるから、本明細書ではこれ以上の詳細な説明は簡略化する。
【００３６】
摩耗進展解析部３２は、車両運動解析部３１により算出された力を用いて模擬レールの摩耗の進展速度を算出する。本実施の形態である摩耗進展解析部３２は、レール接触解析部３２ａ、摩耗量算出部３２ｂ及びレールモデルデータ更新部３２ｃを有する。
【００３７】
レール接触解析部３２ａは、車両運動解析部３１により算出された力を用いて、模擬車両５０と模擬レール５１との接触面内における接触応力を算出する。摩耗量算出部３２ｂは、レール接触解析部３２ａにより算出された接触応力に基づいて、模擬レールの摩耗の進展速度を算出する。レールモデルデータ更新部３２ｃは、摩耗量算出部３２ｂにより算出された摩耗量に基づいてレールモデルデータ４１（レール形状データ４２を含む）を更新する。そして、車両運動解析部３１は、レールモデルデータ更新部３２ｃにより更新されたレールモデルデータ４１を用いて模擬レール５１に作用する力を算出する。
【００３８】
き裂進展解析部３３は、車両運動解析部３１により算出された力を用いて、模擬レール５１に生じうるき裂の進展速度を算出する。本実施の形態であるき裂進展解析部３３は、初期値設定部３３ａ、き裂発生条件判定部３３ｂ及びき裂進展量算出部３３ｃを有する。
【００３９】
初期値設定部３３ａは、後述するき裂発生条件判定部３３ｂによる判定動作に必要なパラメータの初期値を設定する。き裂発生条件判定部３３ｂは、レール接触解析部３２ａにより算出された接触応力に基づいて模擬レール５１にき裂が発生するか否かを判定する。好ましくは、き裂発生条件判定部３３ｂは、初期値設定部３３ａにより設定されたパラメータの初期値を用いて、模擬レール５１にき裂が発生するか否かを判定する。き裂進展量算出部３３ｃは、き裂発生条件判定部３３ｂによりき裂が発生すると判定されたら、模擬レールに生じうるき裂の進展速度を算出する。
【００４０】
きしみ割れ推定部３４は、摩耗進展解析部３２、特に摩耗量算出部３２ｂにより算出された摩耗の進展速度と、き裂進展解析部３３、特にき裂進展量算出部３３ｃにより算出されたき裂の進展速度とを用いて、模擬レールにきしみ割れが生じるか否かを推定する。
【００４１】
制御部２０を構成する各機能部については、後に再度詳細に説明する。
【００４２】
（きしみ割れ推定装置の動作）
次に、図２〜図３のフローチャートを参照して、本実施の形態であるきしみ割れ推定システムＳの動作について説明する。なお、制御部２０を構成する各部の説明について詳述した内容については繰り返しの説明を省略することがある。
【００４３】
図２は、本実施の形態であるきしみ割れ推定装置１０を説明するための図である。きしみ割れ推定装置１０の動作が開始されると、まず、ステップＳ１０において、車両運動解析部３１による、記憶部２１に格納された車両モデルデータ４０及びレールモデルデータ４１を用いた車両運動解析処理、より詳細には、模擬レール５１上を模擬車両が走行した際の模擬レールに作用する力を算出する処理が行われる。
【００４４】
（Ａ）車両運動解析
ステップＳ１０において用いられる模擬車両及び模擬レールの概略を図４に、模擬レールの断面形状、つまりレール形状データ４２の詳細を図５〜図７に示す。模擬車両５０（図４では台車部分のみ図示している）は在来線通勤形車両を模擬し、平均的な走行速度で走行させた。模擬車両５０の車輪形状は修正円弧踏面の設計形状を用い、模擬レール５１の断面形状はＪＩＳ５０ｋｇＮレールのものを採用した。車両運動解析部３１による算出処理では、Ｘ軸を模擬レール５１の延長方向、Ｙ軸を模擬レールのまくらぎ方向、Ｚ軸を鉛直方向として絶対座標系を設定した。
【００４５】
なお、レールの摩耗形状は、軌道諸元のほか、走行する車両条件、とりわけ車輪とレールとが接する面（以下、「車輪踏面」と記す）の影響を大きく受ける。しかし、車輪踏面はレールと同様に走行する距離に応じて摩耗するが、その状況は走行する線区（例えば直線と曲線の割合）によって大きく異なる。本実施の形態である車両運動解析部３１では、摩耗によるレールの形状変化に着目するため、車輪形状については摩耗による形状変化を考慮せず、設計形状の踏面を有する車輪のみが通過するものとした。
【００４６】
模擬レール５１の形状を構成する座標データであるレール形状データ４２は、一例として図５に示すように、レール頭頂面と頭頂面から２５ｍｍ下の側面まで０．４ｍｍ間隔で４００点配置した。レールの摩耗を考慮しない場合は、図６のようにレール形状を一定としてモデル化すればよいが、本実施の形態である車両運動解析部３１のように長手方向にも摩耗が変化する場合は、逐次、断面形状を変化させなければならない。そこで本実施の形態では、レール長手方向にレール断面の着目点を例えば５ｍおきに配置し、隣り合う着目点の断面形状をスプライン補間することで、長手方向のレール断面を生成した（図７）。レール形状データ４２は事前に記憶部２１に格納されていてもよく、ステップＳ１０の処理において作成されてもよい。後述するように、レールモデルデータ更新部３２ｃによりレール形状データ４２の更新が行われる。
【００４７】
（Ｂ）レール接触解析
次に、ステップＳ１１では、摩耗進展解析部３２のレール接触解析部３２ａが、ステップＳ１０において車両運動解析部３１により算出された力を用いて、模擬車両５０と模擬レール５１との接触面内における接触応力を算出する処理を行う。
【００４８】
車両運動解析部３１により、模擬車両５０の車輪と模擬レール５１との接触位置や接触部に作用する力が算出されるが、接触応力の分布については、詳細な結果は得られない。そこで本実施の形態であるレール接触解析部３２ａでは、Ｈｅｒｔｚの接触理論により算出された断面方向の接触半径を用いて、接触面内における応力の分布を算出した。
【００４９】
車両運動解析部３１による解析の結果、接触点Ｃおよび接触点における法線方向の力（以下、「接触力」と記す）Ｎが、図８のように与えられているとする。図８の観測点は、模擬レール５１の断面形状のＹ−Ｚ座標系である。本実施の形態であるレール接触解析部３２ａでは、一例として接触点Ｃを最大１０点まで考慮することとし、それぞれＹ軸の正の方向に向かって１，２，３，…，ｉ，…，１０とする。模擬レール５１の断面形状を構成する観測点を、Ｙ軸の負の方向から１，２，３，…，ｉ，…，４００とする。図８の接触部ＣＰを拡大したものを図９に示す。各接触点Ｃにおいて、接触面と平行な軸をη（ｉ）と定義し、接触力Ｎは接触点Ｃの中心を最大値として楕円形に分布すると考える。
車両運動解析部３１より算出された模擬車両５０の車輪と模擬レール５１の接触力Ｎ（ｉ）を用いて、各接触点の最大応力ＰＭａｘ（ｉ）を表すことが出来る。レールと車輪の接触点Ｃの数をＱと仮定すると、最大応力は以下のように考えられる。
【数１】
2019178944.tif 000003
ここでａ（ｉ）、ｂ（ｉ）は、それぞれ接触部ＣＰの楕円の長径（レール長手方向）、短径（断面方向）の接触半幅である。
【００５０】
接触面に平行な軸η（ｉ）とＹ軸のなす角をθｉとすると、Ｙ軸での接触範囲は以下の式のように表される。
【数２】
2019178944.tif 000004
ここでＣ（ｉ）は車両運動解析部３１により算出される接触点中心である。
【００５１】
接触応力は、接触面内においては図９のように分布することから、次式のように記述できる。
【数３】
2019178944.tif 000005
【００５２】
式（４）により算出された接触応力Ｐ（ｉ，ｊ）は、記憶部２１に応力データ４３として少なくとも一時的に格納される。
【００５３】
（Ｃ）摩耗量算出
次に、ステップＳ１２では、摩耗進展解析部３２の摩耗量算出部３２ｂが、ステップＳ１１においてレール接触解析部３２ａにより算出された接触応力Ｐ（ｉ，ｊ）に基づいて、模擬レール５１の摩耗の進展速度を算出する。
【００５４】
車輪とレールの転がり接触による摩耗は主に凝着摩耗であると考えられ，凝着摩耗を表す摩耗則はこれまでにいくつか提案されている．本実施の形態である摩耗量算出部３２ｂでは，Archardの摩耗予測式(J.F.Archard，“Contact and Rubbing of Flat Surface”，J.Appl.Phys，Vol.24，1953，pp.981-988)を適用した．Archardの摩耗予測式は，次の式（５）で与えられる。
【数４】
2019178944.tif 000006
ただし、Ｗは摩耗体積、Ｆは接触荷重、Ｈは接触する物体のうち柔らかい方の材料硬さ（ここではレール材のビッカーズ硬さ）である。またsはすべり距離である。ｋは摩耗係数であり、材料固有の値である。
【００５５】
式（５）より、模擬レール５１の各接触点Ｃにおける摩耗深さは、次の式（６）で与えられる。
【数５】
2019178944.tif 000007
ここでｄは摩耗深さ、Ｐは接触面圧、δは単位長さあたりのすべり距離（すべり率と等価）であり、本発明者らが行った室内摩耗試験結果より2.54×10-4とした。 
【００５６】
なお、本実施の形態である摩耗量算出部３２ｂにおける摩耗形状予測モデルは、車輪／レール間の巨視的な摩耗に着目するため、接触面内におけるすべり／固着領域の区別は行わず、一定とした。
【００５７】
ステップＳ１２で算出された摩耗量（摩耗深さｄ）は、摩耗量データ４４として記憶部２１に少なくとも一時的に格納される。
【００５８】
次に、ステップＳ１３では、き裂進展解析部３３が、ステップＳ１０において車両運動解析部３１により算出された力を用いて、模擬レール５１に生じうるき裂の進展速度を算出する。
【００５９】
（Ｄ）き裂進展解析
ステップＳ１３におけるき裂進展解析処理の詳細について、図３のフローチャートを参照して説明する。
【００６０】
（Ｄ−１）初期値設定
ステップＳ２０では、き裂進展解析部３３の初期値設定部３３ａが、後述するステップＳ２１において行われる、き裂発生条件判定部３３ｂによる判定動作に必要なパラメータの初期値を設定する。初期値設定部３３ａが設定する初期値とは、模擬レール５１内のフェライト相・パーライト相・非金属材料の分布の計算、結晶粒のゆがみ、き裂の原点、スリップラインの向きに関するパラメータである。初期値設定部３３ａは、これらパラメータを乱数を発生して設定している。
【００６１】
初期値設定部３３ａにより設定されたパラメータの初期値は、き裂発生初期値データ４５として記憶部２１に少なくとも一時的に格納される。
【００６２】
（Ｄ−２）き裂発生条件判定
次に、ステップＳ２１では、き裂進展解析部３３のき裂発生条件判定部３３ｂが、初期値設定部３３ａにより設定されたパラメータの初期値を用いて、ステップＳ１１においてレール接触解析部３２ａにより算出された接触応力Ｐに基づいて模擬レール５１にき裂が発生する条件を算出し、ついでステップＳ２２では、き裂発生条件判定部３３ｂが模擬レール５１にき裂が発生するか否かを判定する。そして、模擬レール５１にき裂が発生したと判定した（ステップＳ２２においてＹＥＳ）場合は、プログラムはステップＳ２３に進み、模擬レール５１にき裂が発生しないと判定したらプログラムはステップＳ１４に戻る。
【００６３】
き裂の種となるフェライト相・パーライト相・非金属材料の分布、及びき裂方向（スリップライン）となる偏向方向（層方向）や境界方向は、ステップＳ２０において初期値設定部３３ａにより乱数で決定されるが、これらのき裂が実際に進展するかどうかは、塑性流動方向に対し定められた閾値角度内にき裂方向が存在するかどうかによる。この条件を満たした場合、その結晶粒のき裂の種は活性化し、き裂発生条件が評価できるようになる。
【００６４】
模擬レール５１におけるレール鋼は、（ａ）フェライト相、（ｂ）パーライト相及び（ｃ）非鉄金属材料から構成される。ここでは各相のき裂の種の活性条件、進展ルールについて説明する。
【００６５】
（ａ）フェライト相
フェライト相は結晶粒内部をき裂が進展する。結晶粒内部にき裂の種を持ち、偏向方向（層方向）をき裂方向（スリップライン）とする。塑性流動方向に対し定められた閾値角度内にき裂方向がある場合、き裂の種は活性化する。
【００６６】
（ｂ）パーライト相
パーライト相は結晶粒内部または境界をき裂が進展する。パーライト相では、き裂の種は境界上に配置する。パーライト相のき裂の種の発生ルールを図示すると図１０のようになる。
【００６７】
パーライト相におけるき裂の種６０は、ランダムに生成した点６１と重心６２を結ぶ直線６３と境界線６４の交点とする（ステップ１）。パーライト相が隣り合う場合において、すでに隣接する要素のき裂の種６０が境界線６４上に存在する場合は、別の境界線６４上にき裂の種６０が生成されるようにする（ステップ２）。
【００６８】
パーライト相のき裂の種の活性と進展は以下のルールに従うものとする。
（１）界面進展き裂：き裂の種が配置された結晶粒界面の方向が、塑性流動方向に対し定められた閾値角度内の場合、き裂の種は活性化し、き裂は結晶粒界面上を進展する。
（２）結晶粒内部進展き裂：結晶粒の偏向方向（層方向）が塑性流動方向に対し定められた閾値角度以内なら、き裂の種は活性し、結晶粒内をき裂が進展する。（１）と（２）とでは（１）が優先される。
（３）パーライト相の結晶粒界面を進展したき裂がフェライト相にぶつかったときは、フェライト相中を進展したき裂と合体するまでは止まっていることとする。
（４）き裂の進展則は、フェライト相と同じものを使う。
【００６９】
（ｃ）非金属材料
非金属材料はパーライト相にのみ存在するものとし、その位置はき裂の種の位置とする。そのため、非金属材料の種の発生条件はパーライト相と同様である。いくつかのパーライト相の中で非金属材料を包含することになったものは、その種が非金属材料の種として利用されるイメージである。
【００７０】
非金属材料のき裂は結晶粒界面を進展する。き裂の種の結晶粒界面の方向が、塑性流動方向に対し定められた閾値角度内の場合、き裂の種は活性化する。
【００７１】
ステップ１におけるき裂の進展ルールはパーライト相と同様であるが、ステップ２の考え方が異なる。累積せん断歪みが大きくなり、変形が進み変形後の非金属の長いほうの対角線が規定値（結晶粒径の３倍）を超えた場合、ステップ２のルールに従いき裂が進展する。
【００７２】
非金属材料におけるステップ２のき裂は、連結されたき裂の先端どうしを結んだ方向をき裂方向とする。この方向が塑性流動方向に対し定められた閾値角度内にない場合は進展しないものとする。
【００７３】
本実施の形態であるき裂発生条件判定部３３ｂでは、非金属材料について累積せん断歪みを算出し、その値から変形量を計算している。図１１に示すように、正方形であったものが、変形によりひし形に歪んだ場合、塑性流動方向に対し定められた閾値角度は見かけ上大きくなる。また、ステップ２に進む際にき裂の先端距離で連結の判定を行うが、変形に伴いこの距離も変わってくる。変形に伴い閾値角度およびき裂先端距離の評価を再度行なうことを、種と距離の動的判定と呼ぶことにする。
【００７４】
閾値角度θは
【数６】
2019178944.tif 000008
で補正させる。先端座標は１次変換行列
【数７】
2019178944.tif 000009
により写像させた点で評価する。
【００７５】
本実施の形態であるき裂発生条件判定部３３ｂにおいて判定されるき裂発生条件は、（ａ）ラチェッティングによるき裂発生条件と、（ｂ）疲労によるき裂発生条件と、（ｃ）これらのうち早いものの３つの条件が存在する。ラチェッティングは累積せん断歪みが規定値を超えた場合の条件であり、変形による発生条件である。疲労は各種モデルで計算された限界サイクル数を実際のサイクル数が超えた場合の条件である。
【００７６】
（ａ）ラチェッティングによるき裂発生条件
（ａ−１）累積せん断ひずみモデル
フェライト、パーライトの各要素には、材料特性（初期せん断降伏応力、破壊の限界せん断ひずみ）を割り当てる。初期せん断降伏応力は以下の式に従う。
[数８]
ｋ０＝0.8×106Ｈｎ                             （７）
ここで、
ｋ０：初期せん断降伏応力
Ｈｎ：ナノ硬さ
である。Ｈｎはフェライト相、パーライト相で以下の値とする。
【表１】
2019178944.tif 000010
【００７７】
破壊の限界せん断ひずみγｃは、フェライト、パーライト共に以下の値とする。
[数９]
γｃ＝11                                       （８）
【００７８】
本実施の形態であるき裂発生条件判定部３３ｂでは、荷重サイクルが進むにつれてせん断ひずみが累積するモデルを採用する。モデルは、以下のものを採用する。
【数１０】
2019178944.tif 000011
ここで、
【数１１】
2019178944.tif 000012
である。
【００７９】
式（７）により求められたフェライト・パーライトの各要素における全累積せん断ひずみγｉｊが、要素毎に初期に割り当てられた破壊の限界せん断ひずみγｃを超えた場合に、き裂が発生するものとする。
【００８０】
また、パラメータα、βはフェライト、パーライトで以下の値とする。
【表２】
2019178944.tif 000013
【００８１】
累積せん断ひずみは、レール接触解析部３２ａにより算出された車輪／レール接触の応力値より算出する。累積せん断ひずみの算出には最大直交せん断応力τｊｚｘ（ｍａｘ）が必要であるが、この応力値はレール接触解析部３２ａによる車輪／レール接触解析の結果を用いる。
【００８２】
ここで、累積せん断ひずみは粒界フェライト部分についても考慮している。フェライト相の場合、その結晶粒の累積せん断ひずみはフェライトの材料特性より算出したものを用いるが、パーライト相の場合は、粒界と結晶粒のそれぞれの材料特性より求めた累積せん断ひずみに体積比を乗じて合計する。
【００８３】
ここで、パーライトの累積せん断歪みをrp、フェライト相の累積せん断歪みをrfとすると、全体の累積せん断歪みは
[数１２]
rp×(1-Ab/Aa)+rf×Ab/Aa                       （１２）
となる。このような計算で求めた累積せん断歪みを各結晶粒に保持する。ラチェッティングによるき裂の発生判定には、上記で求めた全体の累積せん断歪みの値を利用する。
【００８４】
（ａ−２）累積垂直ひずみモデル
き裂の発生条件には関係しないが、垂直方向の変形算出に際して累積せん断ひずみと同様の方法を用いている。
【数１３】
2019178944.tif 000014
上記の式は式（７）〜（１１）に対応したものである。これらの式からある行（同じ高さの粒子群）のひずみを求めることができ、車輪通過ごとの垂直方向の変位は、Δεｉｊ×（ある行の高さ）で求められる。
【００８５】
（ｂ）疲労によるき裂発生条件
疲労によるき裂発生は、荷重サイクル数が下式（１７）で得られるＮを超えた際に発生するものとする。
【数１４】
2019178944.tif 000015
ここで、すべり帯の長さｄは、図１２に示すように、各結晶粒において、重心６２を通り塑性流動方向に伸ばした直線６３がそれぞれ界面６４と交わる線分の長さとする。
【００８６】
また、下式（１８）により、分解せん断応力範囲Δτｒｅｓを算出する。
[数１５]
Δτｒｅｓ＝τｍａｘ−τｍｉｎ                        （１８）
ここで、τｍａｘ、τｍｉｎはレール接触解析部３２ａによる車輪／レール接触解析の結果から取得することができる。応力の評価の座標は仮想的なき裂の種、方向は仮想的なき裂進展方法とする。
【００８７】
（Ｄ−３）き裂進展量算出
ステップＳ２３では、き裂発生条件判定部３３ｂによりき裂が発生すると判定されたので（ステップＳ２２においてＹＥＳ）、き裂進展解析部３３のき裂進展量算出部３３ｃが模擬レール５１に生じうるき裂の進展速度を算出する。
【００８８】
き裂進展量算出部３３ｃにおいて算出されるき裂進展量におけるき裂進展規則は、フェライト相・パーライト相・非金属材料に共通である。
【００８９】
本実施の形態であるき裂進展量算出部３３ｃにおけるき裂進展のアルゴリズムは、H. A. Suharutono，K. Potter，A. Schram and H. Zenner，“Multiaxial Fatigue and Deformation: Testing and Prediction”，ASTM STP 1387，2000，pp.323を基にしている。
【００９０】
（ａ）第1ステップ
結晶粒内を進展するき裂は、図１３に示すように、き裂６５の原点（き裂の種）６０から、スリップライン（き裂方向）６６に沿って式（１９）に従って進展する。ここで、ａはき裂の長さであり、Ｎは応力の負荷回数である。
【数１６】
2019178944.tif 000016
式（２０）は、スリップラインに沿って働くせん断応力の変動幅を表しており、Δσｘ、Δσｙ、ΔσｘｙはＸ、Ｙ軸で定義される垂直応力とせん断応力の変動幅である。き裂進展量算出部３３ｃでは、き裂発生条件判定部３３ｂで求めた値を用いてき裂進展量を算出する。
【００９１】
（ｂ）第２ステップ
複数のき裂は、その先端の距離が結晶粒径の２５％以内であり、且つ、き裂自身の長さが結晶粒径の７５％に達すると互いに結合する。ただし、結合直前のき裂の形状は、結合後に変化することはない。また、き裂の長さが結晶粒径の３倍に達したとき今度は、き裂は次式（２１）に従って進展する。
【数１７】
2019178944.tif 000017
ただし、複数のき裂が結合した場合、き裂の長さを、最も長くなる先端同士の直線距離で定義する。また、き裂面はその直線方向で定義し、θはその角度である。き裂の角度が塑性流動方向に対し定められた閾値角度内に含まれない場合は進展しない。Δεは、この面にかかる垂直歪みの変動幅である。
【００９２】
この段階では、き裂はき裂面に対する垂直歪みが引張時においてのみ進展するため、式（２２）の歪みの値が負になった場合は、計算上Δεは０で置き換えられる。また、式（１９）で計算されるき裂の進展が式（２１）の結果よりも大きい場合には、式（１９）の値を用いている。
【００９３】
本実施の形態であるき裂進展量算出部３３ｃでは、初期構造として正六角形状からの結晶粒のゆがみを与えることができる。このゆがみ形状は、星出敏彦、河端享介、井上達雄、日本機械学会論文集Ａ編、５５巻、１９８９年、pp.222を基にしたアルゴリズムを用いて与えられ、図１４に示すように、六角形の各頂点を０＜ｒ＜ζＬの範囲内でランダムに移動させている。ここでＬは正六角形の一辺の長さであり、各結晶粒の大きさを決める。ζは、歪みの大きさを表すパラメータである。
【００９４】
結晶粒を正六角形から歪ませると、結晶粒径を正しく定義することはできない。そのため、き裂の結合判定やき裂の進展速度式の選択に使われる結晶粒径は、本実施の形態では、結晶粒を歪ませる前の半径の２倍(つまり２Ｌ)として定義している。
【００９５】
本実施の形態であるき裂進展量算出部３３ｃでは、以下の摩耗モデルを実装している。摩耗された領域においてき裂が進展することは物理的にありえないので、き裂進展量算出部３３ｃにおいて、摩耗された領域のき裂が進展しない処理が行われている。
【００９６】
上述したArchard摩耗モデルの摩耗深さは以下の式で求めることができる。
【数１８】
2019178944.tif 000018
この摩耗深さΔｚ（ｘ，ｙ）は接触面が一要素分移動する際に摩耗する量である。車輪の接触面は、ある面積を持って移動するため、車輪通過１サイクルにおいて一要素が摩耗する量は、接触面積分を重ね合わせた量となる。
【００９７】
また、本実施の形態であるき裂進展量算出部３３ｃでは、粗さ接触モデルを実装している。粗さ接触は、以下の式を直交せん断応力(垂直応力)に乗じることによって定義する（第１ステップおよび第２ステップ共通）。
【数１９】
2019178944.tif 000019
ここで、ｚは接触面からの深さ、Ａ＝２、ｄ＝１５μｍ（影響深さ）とする。
【００９８】
ステップＳ１３及び図３に示すステップにおいて算出されたき裂進展量は、き裂進展量データ４６として記憶部２１に少なくとも一時的に格納される。
【００９９】
（Ｅ）きしみ割れ推定
図２に戻って、ステップＳ１４では、きしみ割れ推定部３４が、摩耗進展解析部３２、特に摩耗量算出部３２ｂにより算出された摩耗の進展速度と、き裂進展解析部３３、特にき裂進展量算出部３３ｃにより算出されたき裂の進展速度とを用いて、模擬レール５１にきしみ割れが生じるか否かを推定する。そして、模擬レール５１にきしみ割れが生じると推定したら（ステップＳ１４においてＹＥＳ）、ステップＳ１５においてきしみ割れ推定部３４が模擬レール５１にきしみ割れが推定する旨の表示を表示装置１２に表示する。一方、模擬レール５１にきしみ割れが生じないと推定したら（ステップＳ１４においてＮＯ）、プログラムはステップＳ１６に進む。
【０１００】
（Ｆ）レールモデルデータ更新
ステップＳ１６において、摩耗進展解析部３２のレールモデルデータ更新部３２ｃは、摩耗量算出部３２ｂにより算出された摩耗量に基づいてレールモデルデータ４１（含むレール形状データ４２）を更新する。そして、プログラムはステップＳ１０に戻り、車両運動解析部３１は、レールモデルデータ更新部３２ｃにより更新されたレールモデルデータ４１を用いて模擬レール５１に作用する力を算出する。
【０１０１】
（本実施の形態の効果）
このように構成された本実施の形態のきしみ割れ推定装置１０では、車両モデルデータ４０及びレールモデルデータ４１を用いて模擬レール５１上を模擬車両５０が走行した際の模擬レール５１に作用する力を算出する車両運動解析部３１と、車両運動解析部３１により算出された力を用いて模擬レール５１の摩耗の進展速度を算出する摩耗進展解析部３２と、車両運動解析部３１により算出された力を用いて模擬レール５１に生じうるき裂の進展速度を算出するき裂進展解析部３３と、摩耗の進展速度とき裂の進展速度とを用いて模擬レール５１にきしみ割れが生じるか否かを推定するきしみ割れ推定部３４とを有する。
【０１０２】
このようにすることで、模擬レール５１と模擬車両５０の車輪との接触状態に基づいてきしみ割れの発生を定量的に推定することができる。
【０１０３】
ここで、摩耗進展解析部３２が、車両運動解析部３１により算出された力を用いて模擬車両５０と模擬レール５１との接触面内における接触応力を算出するレール接触解析部３２ａと、レール接触解析部３２ａにより算出された接触応力に基づいて模擬レール５１の摩耗の進展速度を算出する摩耗量算出部３２ｂとを有するので、摩耗進展解析部３２による模擬レール５１の摩耗の進展速度算出処理をより正確に定量的に行うことができ、結果として、きしみ割れの発生をより定量的に推定することができる。
【０１０４】
また、摩耗進展解析部３２が摩耗量算出部３２ｂにより算出された摩耗量に基づいてレールモデルデータ４１を更新するレールモデルデータ更新部３２ｃを有し、車両運動解析部３１がレールモデルデータ更新部３２ｃにより更新されたレールモデルデータ４１を用いて模擬レール５１に作用する力を算出するので、模擬レール５１の摩耗量を的確にフィードバックして模擬レール５１に作用する力を算出することができ、結果として、きしみ割れの発生をより定量的に推定することができる。
【０１０５】
さらに、き裂進展解析部３３が、レール接触解析部３２ａにより算出された接触応力に基づいて模擬レール５１にき裂が発生するか否かを判定するき裂発生条件判定部３３ｂと、き裂発生条件判定部３３ｂによりき裂が発生すると判定されたら模擬レール５１に生じうるき裂の進展速度を算出するき裂進展量算出部３３ｃとを有するので、き裂進展解析部３３による、模擬レール５１に生じうるき裂の進展速度算出処理をより正確に定量的に行うことができ、結果として、きしみ割れの発生をより定量的に推定することができる。
【０１０６】
そして、き裂進展解析部３３がき裂発生条件判定部３３ｂによる判定動作に必要なパラメータの初期値を設定する初期値設定部３３ａを有し、き裂発生条件判定部３３ｂが初期値設定部３３ａにより設定されたパラメータの初期値を用いて模擬レール５１にき裂が発生するか否かを判定するので、模擬レール５１と模擬車両５０の車輪との接触状態に基づいてきしみ割れの発生を定量的に推定することができる。
【０１０７】
以上、図面を参照して、本発明の実施の形態を詳述してきたが、具体的な構成は、この実施の形態及び実施例に限らず、本発明の要旨を逸脱しない程度の設計的変更は、本発明に含まれる。
【０１０８】
例えば、上述の実施の形態であるきしみ割れ推定システムＳでは、推定装置１０内に記憶部２１を設けていたが、記憶部２１を推定装置１０とは別体に構成することも可能である。
【０１０９】
そして、上述の実施例において、推定装置１０を動作させるプログラムは記憶部２１に格納されて提供されていたが、不図示の光学ディスクドライブ等を用いて、プログラムが格納されたＤＶＤ（Digital Versatile Disc）、ＵＳＢ外部記憶装置、メモリーカード等を接続し、このＤＶＤ等からプログラムを推定装置１０に読み込んで動作させてもよい。また、インターネット上のサーバ装置内にプログラムを格納しておき、推定装置１０に通信部を設けてこのプログラムを推定装置１０に読み込んで動作させてもよい。さらに、上述の実施例において、推定装置１０は複数のハードウェア要素により構成されていたが、これらハードウェア要素の一部の動作を制御部２０がプログラムの動作により実現することも可能である。
【符号の説明】
【０１１０】
Ｓ    推定システム
１０    推定装置
２０    制御部
２１    記憶部
３０    表示制御部
３１    車両運動解析部
３２    摩耗進展解析部
３２ａ  レール接触解析部
３２ｂ  摩耗量算出部
３２ｃ  レールモデルデータ更新部
３３    き裂進展解析部
３３ａ  初期値設定部
３３ｂ  き裂発生条件判定部
３３ｃ  き裂進展量算出部
３４    きしみ割れ推定部
４０    車両モデルデータ
４１    レールモデルデータ
５０    模擬車両
５１    模擬レール
【図１】
2019178944.tif 000020
【図２】
2019178944.tif 000021
【図３】
2019178944.tif 000022
【図４】
2019178944.tif 000023
【図５】
2019178944.tif 000024
【図６】
2019178944.tif 000025
【図７】
2019178944.tif 000026
【図８】
2019178944.tif 000027
【図９】
2019178944.tif 000028
【図１０】
2019178944.tif 000029
【図１１】
2019178944.tif 000030
【図１２】
2019178944.tif 000031
【図１３】
2019178944.tif 000032
【図１４】
2019178944.tif 000033