(57)【要約】【課題】弾性状態から塑性状態に遷移する構造物に対しても正確な解析を可能とする構造解析装置等を提供する。【解決手段】構造解析装置１は、モデルデータ２の各粒子に対して加速度を計算し、この加速度から時間増分経過後の歪増分を計算する材料状態計算手段３と、記憶された各粒子の状態と歪増分とから、時間増分経過後の当該粒子が弾性状態であるか塑性状態であるかを判断する材料状態判定手段４と、当該粒子の状態が弾性状態である場合にその粒子の弾性応力を計算して降伏したか否かを判断する降伏判定手段５と、粒子が降伏した場合に当該粒子の応力を修正計算して当該粒子の状態を修正記憶する応力修正計算手段６と、記憶された各粒子の状態を外部に出力する終了判定処理手段７とを有する。【選択図】図１
2014137747.tif 000002
【特許請求の範囲】
【請求項１】
構造物を複数の粒子でモデル化し、前記各粒子の状態を含むモデルデータを予め記憶し、このモデルデータに対して粒子法による構造解析を行う構造解析装置であって、
前記モデルデータの前記各粒子に対して加速度を計算し、この加速度から時間増分経過後の歪増分を計算する材料状態計算手段と、
記憶された前記各粒子の状態と前記歪増分とから、前記時間増分経過後の当該粒子が弾性状態であるか塑性状態であるかを判断する材料状態判定手段と、
当該粒子の状態が弾性状態である場合にその粒子の弾性応力を計算して降伏したか否かを判断する降伏判定手段と、
前記粒子が降伏した場合に当該粒子の応力を修正計算して当該粒子の状態を修正記憶する応力修正計算手段と、
記憶された前記各粒子の状態を外部に出力する終了判定処理手段と
を有することを特徴とする構造解析装置。
【請求項２】
前記終了判定処理手段が、
対象となる全ての粒子および時間に対して処理が終了したか否かを判断し、この処理が終了したと判断された場合に記憶された前記各粒子の状態を外部に出力することを特徴とする、請求項１に記載の構造解析装置。
【請求項３】
前記降伏判定手段が、当該粒子の試行相当応力を計算することによってその粒子が降伏したか否かを判断する機能を有すると共に、
前記応力修正計算手段が、前記粒子が降伏した場合に前記試行相当応力に基づいて当該粒子の応力を修正計算する機能を有することを特徴とする、請求項１に記載の構造解析装置。
【請求項４】
前記材料状態計算手段が、前記各粒子に対して近傍粒子を探索して異材料間の接触の有無を判定し、接触がある場合に離散化された運動方程式から算出された加速度を当該粒子の加速度に加算する機能を有することを特徴とする、請求項１に記載の構造解析装置。
【請求項５】
前記降伏判定手段が、当該粒子の状態が塑性状態である場合に、その粒子の弾塑性応力を計算して降伏したか否かを判断する機能を有することを特徴とする、請求項１に記載の構造解析装置。
【請求項６】
構造物を複数の粒子でモデル化し、前記各粒子の状態を含むモデルデータを予め記憶し、このモデルデータに対して粒子法による構造解析を行う構造解析装置にあって、
前記モデルデータの前記各粒子に対して材料状態計算手段が加速度を計算し、
この加速度から時間増分経過後の歪増分を前記材料状態計算手段が計算し、
記憶された前記各粒子の状態に材料状態判定手段が前記歪増分を加算し、
前記時間増分経過後の当該粒子が弾性状態であるか塑性状態であるかを前記材料状態判定手段が判断し、
当該粒子の状態が弾性状態である場合に降伏判定手段がその粒子の弾性応力を計算し、
この弾性応力からその粒子が降伏したか否かを前記降伏判定手段が判断し、
前記粒子が降伏した場合に応力修正計算手段が当該粒子の応力を修正計算して当該粒子の状態を修正記憶し、
記憶された前記各粒子の状態を終了判定処理手段が外部に出力する
ことを特徴とする構造解析方法。
【請求項７】
構造物を複数の粒子でモデル化し、前記各粒子の状態を含むモデルデータを予め記憶し、このモデルデータに対して粒子法による構造解析を行う構造解析装置にあって、
前記構造解析装置が備えるプロセッサに、
前記モデルデータの前記各粒子に対して加速度を計算する手順、
この加速度から時間増分経過後の歪増分を計算する手順、
記憶された前記各粒子の状態に前記歪増分を加算する手順、
前記時間増分経過後の当該粒子が弾性状態であるか塑性状態であるかを判断する手順、
当該粒子の状態が弾性状態である場合にその粒子の弾性応力を計算する手順、
この弾性応力からその粒子が降伏したか否かを判断する手順、
前記粒子が降伏した場合に当該粒子の応力を修正計算して当該粒子の状態を修正記憶する手順、
および記憶された前記各粒子の状態を終了判定処理手段が外部に出力する手順
を実行させることを特徴とする構造解析プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、構造解析装置、構造解析方法、および構造解析プログラムに関し、特に構造物が弾性状態から塑性状態に遷移する際にも正確な解析を可能とする構造解析装置等に関する。
【背景技術】
【０００２】
構造物の変形、たとえば金属の切削加工やプレス加工などにおいて加工が行われている最中の当該構造物の挙動をシミュレーションする手法の一つとして、有限要素法が知られている。その一例として、非特許文献１には金属板の切削加工に対して有限要素法を適用し弾塑性解析を行うための手法が記載されている。
【０００３】
金属は、切削加工が進むと変形して、弾性状態から永久変形が残る塑性状態に遷移する。その際、弾性の応力増分−歪増分の関係式から応力増分を算出してこれを時間増分前の応力に加算し、現時間の応力を算出する。弾塑性理論では応力は降伏曲面の方程式を満たすと考えられているが、弾性の応力増分−歪増分の関係式から算出した応力はその方程式を満たさない。そのため、状態の遷移に応じて歪や応力を修正する必要がある。
【０００４】
一方、そのような構造物の挙動をシミュレーションする手法で、近年活発に研究がなされ注目されている手法の一つとして、粒子法がある。その中でも特に、ＳＰＨ（Smoothed Particle Hydrodynemics）やＭＰＳ（Moving Particle Simulation）などの手法が広く用いられている。
【０００５】
粒子法は、大きさが同一の複数の粒子によって構造物をモデル化するので、要素単位でモデルを作成しなければならない有限要素法に比べて、複雑な形状の構造物でも短時間で容易にモデルを作成できるという利点がある。また、構造物が大きく変形する場合には、有限要素法では要素が潰れるので解析のために要素作成を何回もやり直すことが必要になるが、粒子法ではその必要がないという利点もある。
【０００６】
これに関連する技術文献として、以下の各々がある。その中でも特許文献１には、シェル要素モデルによって金属のプレス加工をシミュレートするという解析方法が記載されている。特許文献２には、第１及び第２構造体の接触による変形を、粒子モデルによってシミュレートするという方法が記載されている。特許文献３には、アルミニウムの鋳塊形状の予測を凝固温度解析と弾塑性応力解析によって行うという技術が記載されている。
【０００７】
特許文献４には、動的陽解法有限要素法で応力の計算結果の過度な振動を抑制するという技術が記載されている。特許文献５には、多孔質部材が破壊応力に達した場合に要素応力＝０とするという非線形計算解析方法が記載されている。特許文献６には、弾性応力テンソルから弾性ストレス状態を算出するという材料解析方法が記載されている。
【０００８】
特許文献７には、非線形固体問題を有する対象物に対しても重合メッシュ法を適用可能という有限要素解析方法が記載されている。特許文献８には、人体の血管内部の血流を分析するという超音波診断装置が記載されている。特許文献９には、構造物を複数の要素に分割して有限要素法によって解析するという構造解析装置が記載されている。
【先行技術文献】
【特許文献】
【０００９】
【特許文献１】特開２００４−０４２０９８号公報
【特許文献２】特開２００８−１５２４２３号公報
【特許文献３】特開平１１−０５７９４７号公報
【特許文献４】特開２００１−２０１４０８号公報
【特許文献５】特開２００５−３０８５９５号公報
【特許文献６】特開２００６−３１３０９４号公報
【特許文献７】特開２００９−０８６８０７号公報
【特許文献８】特開２０１０−１２５２０３号公報
【特許文献９】特開２０１０−２３０６４４号公報
【非特許文献】
【００１０】
【非特許文献１】真鍋他、「陰解法弾塑性有限要素法による切削機構の解析−radial return法による定式化と切りくずカール解析−」、精密工学会誌、ｖｏｌ．６２、Ｎｏ．８、１９９６年
【発明の概要】
【発明が解決しようとする課題】
【００１１】
金属の切削加工など、構造物の弾塑性特性を考慮する必要がある構造解析においても、粒子法を適用できれば、上記のような利点を生かして高精度のシミュレーションを短時間で行うことができるであろうと考えられており、実際にそのような試みも一部で行われている。
【００１２】
しかしながら、既存技術に係る粒子法では、構造物の弾性状態から塑性状態への遷移は考慮されていない。構造物の変形に伴ってこの遷移が行われる場合に、前述の非特許文献１に記載の有限要素法と同様に、歪や応力の修正が必要となる。この修正がないと解析精度が大きく低下することとなる。
【００１３】
この点を解決し、弾性状態から塑性状態への遷移が行われる構造物に対しても粒子法を適用した構造解析を可能とする技術は、前述の特許文献１〜９および非特許文献１には記載されていない。
【００１４】
本発明の目的は、弾性状態から塑性状態に遷移する構造物に対しても粒子法を適用して正確な解析を可能とする構造解析装置、構造解析方法、および構造解析プログラムを提供することにある。
【課題を解決するための手段】
【００１５】
上記目的を達成するため、本発明に係る構造解析装置は、構造物を複数の粒子でモデル化し、各粒子の状態を含むモデルデータを予め記憶し、このモデルデータに対して粒子法による構造解析を行う構造解析装置であって、モデルデータの各粒子に対して加速度を計算し、この加速度から時間増分経過後の歪増分を計算する材料状態計算手段と、記憶された各粒子の状態と歪増分とから、時間増分経過後の当該粒子が弾性状態であるか塑性状態であるかを判断する材料状態判定手段と、当該粒子の状態が弾性状態である場合にその粒子の弾性応力を計算して降伏したか否かを判断する降伏判定手段と、粒子が降伏した場合に当該粒子の応力を修正計算して当該粒子の状態を修正記憶する応力修正計算手段と、記憶された各粒子の状態を外部に出力する終了判定処理手段とを有することを特徴とする。
【００１６】
上記目的を達成するため、本発明に係る構造解析方法は、構造物を複数の粒子でモデル化し、各粒子の状態を含むモデルデータを予め記憶し、このモデルデータに対して粒子法による構造解析を行う構造解析装置にあって、モデルデータの各粒子に対して材料状態計算手段が加速度を計算し、この加速度から時間増分経過後の歪増分を材料状態計算手段が計算し、記憶された各粒子の状態に材料状態判定手段が歪増分を加算し、時間増分経過後の当該粒子が弾性状態であるか塑性状態であるかを材料状態判定手段が判断し、当該粒子の状態が弾性状態である場合に降伏判定手段がその粒子の弾性応力を計算し、この弾性応力からその粒子が降伏したか否かを降伏判定手段が判断し、粒子が降伏した場合に応力修正計算手段が当該粒子の応力を修正計算して当該粒子の状態を修正記憶し、記憶された各粒子の状態を終了判定処理手段が外部に出力することを特徴とする。
【００１７】
上記目的を達成するため、本発明に係る構造解析プログラムは、構造物を複数の粒子でモデル化し、各粒子の状態を含むモデルデータを予め記憶し、このモデルデータに対して粒子法による構造解析を行う構造解析装置にあって、構造解析装置が備えるプロセッサに、モデルデータの各粒子に対して加速度を計算する手順、この加速度から時間増分経過後の歪増分を計算する手順、記憶された各粒子の状態に歪増分を加算する手順、時間増分経過後の当該粒子が弾性状態であるか塑性状態であるかを判断する手順、当該粒子の状態が弾性状態である場合にその粒子の弾性応力を計算する手順、この弾性応力からその粒子が降伏したか否かを判断する手順、粒子が降伏した場合に当該粒子の応力を修正計算して当該粒子の状態を修正記憶する手順、および記憶された各粒子の状態を終了判定処理手段が外部に出力する手順を実行させることを特徴とする。
【発明の効果】
【００１８】
本発明は、上述したように、当該粒子が降伏したか否かを判断して、降伏した場合にその粒子の応力を修正計算するよう構成したので、降伏による弾性状態から塑性状態への遷移を含めた解析が可能となる。これによって、弾性状態から塑性状態に遷移する構造物に対しても粒子法を適用して正確な解析が可能であるという、優れた特徴を持つ構造解析装置、構造解析方法、および構造解析プログラムを提供することができる。
【図面の簡単な説明】
【００１９】
【図１】本発明の基本形態に係る構造解析装置の構成について示す説明図である。
【図２】本発明の実施形態に係る構造解析装置の構成について示す説明図である。
【図３】図２に示した構造解析装置の動作について示すフローチャートである。
【図４】図３の続きである。
【図５】図３のステップＳ１１２〜１１３として示した降伏判定手段および応力修正計算手段による処理の詳細について示すフローチャートである。
【図６】図２〜４に示した構造解析装置で行う計算の対象となる粒子モデルについて示す説明図である。
【図７】図６に示した粒子モデルで、板に衝突する直方体の接触面の中央部に位置する粒子に関して、相当塑性歪と相当応力の関係を図２〜４に示した構造解析装置でシミュレートした結果について示すグラフである。
【発明を実施するための形態】
【００２０】
（基本形態）
以下、本発明の基本形態の構成について添付図１に基づいて説明する。
最初に、基本形態の基本的な内容について説明し、その後でより具体的な内容について説明する。
基本形態に係る構造解析装置１は、構造物を複数の粒子でモデル化し、各粒子の状態を含むモデルデータ２を予め記憶し、このモデルデータに対して粒子法による構造解析を行う構造解析装置であり、モデルデータの各粒子に対して加速度を計算し、この加速度から時間増分経過後の歪増分を計算する材料状態計算手段３と、記憶された各粒子の状態と歪増分とから、時間増分経過後の当該粒子が弾性状態であるか塑性状態であるかを判断する材料状態判定手段４と、当該粒子の状態が弾性状態である場合にその粒子の弾性応力を計算して降伏したか否かを判断する降伏判定手段５と、粒子が降伏した場合に当該粒子の応力を修正計算して当該粒子の状態を修正記憶する応力修正計算手段６と、記憶された各粒子の状態を外部に出力する終了判定処理手段７とを有する。
【００２１】
以上の構成を備える事により、この構造解析装置１は、構造物が弾性状態から塑性状態に遷移する際にも正確な解析を行うことが可能となる。
以下、これをより詳細に説明する。
【００２２】
図１は、本発明の基本形態に係る構造解析装置１の構成について示す説明図である。構造解析装置１には、構造物を複数の粒子でモデル化して各粒子の状態を含むモデルデータ２が予め記憶されている。そして、このモデルデータ２に記憶された各粒子の状態に対して、材料状態計算手段３、材料状態判定手段４、降伏判定手段５、応力修正計算手段６、終了判定処理手段７が順次計算を行い、その結果を外部に対して出力する。
【００２３】
材料状態計算手段３は、各粒子に対して加速度を計算し、この加速度から時間増分経過後の歪増分を計算する。材料状態判定手段４は、記憶された各粒子の状態と歪増分とから、時間増分経過後の当該粒子が弾性状態であるか塑性状態であるかを判断する。その粒子の状態が弾性状態である場合に、降伏判定手段５はその粒子の弾性応力を計算して降伏したか否かを判断する。
【００２４】
そして、粒子が降伏した場合に、応力修正計算手段６が当該粒子の応力を修正計算して当該粒子の状態を修正記憶する。最終的に記憶された各粒子の状態を、終了判定処理手段７が外部に対して出力する。
【００２５】
この構成を備えることによって、弾性状態から塑性状態に遷移する際の応力が反映された粒子の状態を出力し、より正確な解析結果を得ることが可能となる。なお、上記各処理の具体的な計算内容については、この後の実施形態で詳しく述べる。
【００２６】
（実施形態）
続いて、本発明の実施形態の構成について添付図２に基づいて説明する。
最初に、本実施形態の基本的な内容について説明し、その後でより具体的な内容について説明する。
本実施形態に係る構造解析装置１０は、構造物を複数の粒子でモデル化し、各粒子の状態を含むモデルデータ２７を予め記憶し、このモデルデータに対して粒子法による構造解析を行う構造解析装置である。この構造解析装置１０は、モデルデータの各粒子に対して加速度を計算し、この加速度から時間増分経過後の歪増分を計算する材料状態計算手段２１と、記憶された各粒子の状態と歪増分とから、時間増分経過後の当該粒子が弾性状態であるか塑性状態であるかを判断する材料状態判定手段２２と、当該粒子の状態が弾性状態である場合にその粒子の弾性応力を計算して降伏したか否かを判断する降伏判定手段２３と、粒子が降伏した場合に当該粒子の応力を修正計算して当該粒子の状態を修正記憶する応力修正計算手段２４と、記憶された各粒子の状態を外部に出力する終了判定処理手段２５とを有する。
【００２７】
また、終了判定処理手段２５は、対象となる全ての粒子および時間に対して処理が終了したか否かを判断し、この処理が終了したと判断された場合に記憶された各粒子の状態を外部に出力する。そして降伏判定手段２３は、当該粒子の試行相当応力を計算することによってその粒子が降伏したか否かを判断する機能を有し、応力修正計算手段２４が、粒子が降伏した場合に試行相当応力に基づいて当該粒子の応力を修正計算する機能を有する。
【００２８】
さらに材料状態計算手段２１は、各粒子に対して近傍粒子を探索して異材料間の接触の有無を判定し、接触がある場合に離散化された運動方程式から算出された加速度を当該粒子の加速度に加算する機能を有する。そして降伏判定手段２３は、当該粒子の状態が塑性状態である場合に、その粒子の弾塑性応力を計算して降伏したか否かを判断する機能を有する。
【００２９】
以上の構成を備える事により、この構造解析装置１は、構造物が弾性状態から塑性状態に遷移する際にも、粒子法を適用して正確な解析を行うことが可能となる。
以下、これをより詳細に説明する。
【００３０】
図２は、本発明の実施形態に係る構造解析装置１０の構成について示す説明図である。構造解析装置１０は、一般的なコンピュータ装置としての構成を備えるものである。即ち、コンピュータプログラムを実行する主体であるプロセッサ１１と、処理中のデータを一時的に記憶する主記憶手段１２（ＲＡＭ：Random Access Memory）と、プログラムやデータを非一時的に記憶する固定記憶手段１３とを備える。固定記憶手段１３は、典型的にはハードディスクドライブやフラッシュメモリなどである。
【００３１】
構造解析装置１０はさらに、処理対象のデータを外部から入力する入力手段１４と、処理結果を外部に出力する出力手段１５とを備える。典型的には、入力手段１４はキーボードやマウスなどであり、出力手段１５はディスプレイやプリンタなどであるが、ネットワークを介して接続された他の装置を入力手段１４および出力手段１５として機能させることももちろん可能である。
【００３２】
プロセッサ１１は、塑性変形解析プログラムが動作することにより、後述する材料状態計算手段２１、材料状態判定手段２２、降伏判定手段２３、応力修正計算手段２４、終了判定処理手段２５の各々として機能する。また、主記憶手段１２には後述する材料状態識別情報２６が一時的に記憶され、固定記憶手段１３には解析に使用されるモデルデータ２７が予め記憶されている。モデルデータ２７は、より具体的には粒子の粒子番号、粒子の座標、粒子の速度、粒子の加速度、粒子の材料番号、材料特性値などのような各々の数値を含む。また、モデルデータ２７には処理後の材料状態識別情報２６が反映される。
【００３３】
（解析の基本概念）
本実施形態では、粒子法の中でもＳＰＨ（Smoothed Particle Hydrodynamics）を用いた解析を実施する。公知のように、ＳＰＨは構造物を複数の粒子でモデル化し、その中でも注目する粒子の物理量については重み関数を用いて近似的に表現する。以下、その基本概念について説明する。
【００３４】
重み関数として、３次スプライン関数を用いることができる。ｘを注目する粒子の中心の位置、ｘ’を注目する粒子の近傍粒子の中心の位置、ｈを粒子の直径、Ｄを空間次元、さらにｚを以下の数１のように定義する。
【００３５】
【数１】
2014137747.tif 000003
【００３６】
重み関数Ｗ（ｚ）は、以下の数２のように定義することができる。
【００３７】
【数２】
2014137747.tif 000004
【００３８】
注目する粒子の物理量ｆ（ｘ）は、この重み関数Ｗを用いて、以下の数３のように近似的に表現することができる。
【００３９】
【数３】
2014137747.tif 000005
【００４０】
さらに、Ｊを近傍粒子の番号、ｘＪを粒子Ｊの中心の位置、ｍＪを粒子Ｊの質量、ρＪを粒子Ｊの密度とすると、注目する粒子の物理量ｆ（ｘ）は、以下の数４のように離散化して表現することができる。
【００４１】
【数４】
2014137747.tif 000006
【００４２】
また、ｆ（ｘ）の勾配および発散は各々、以下の数５および数６のように表現することができる。
【００４３】
【数５】
2014137747.tif 000007
【００４４】
【数６】
2014137747.tif 000008
【００４５】
（フローチャート）
図３〜４は、図２に示した構造解析装置１０の動作について示すフローチャートである。図３〜４は、紙面の都合で、１枚のフローチャートを２枚に分けて示している。まず、材料状態計算手段２１がモデルデータ２７を読み込み、この中から必要な情報を主記憶手段１２に一時保存する（ステップＳ１０１）。この中には、後述の材料状態識別情報２６も含まれる。
【００４６】
次に、材料状態計算手段２１は、モデルデータ２７の中のｉ番目の粒子である粒子Ｉを中心とする半径の球内部に位置する粒子を探索し、粒子Ｉの粒子番号と、球内部に位置する粒子Ｉ以外の粒子の粒子番号を関連付けて主記憶手段１２に一時保存する（ステップＳ１０２）。
【００４７】
続いて材料状態計算手段２１は、粒子Ｉの加速度を計算して、粒子Ｉの粒子番号とその加速度とを関連づけて主記憶手段１２に一時保存する（ステップＳ１０３）。そして材料状態計算手段２１は、異材料間の接触の有無を判定する（ステップＳ１０４）。より具体的には、粒子Ｉの中心と異材料の粒子Ｊの中心の距離が粒子の直径以下であればその２粒子は接触したと判定する。接触したらステップＳ１０５に、接触していなければステップＳ１０６に各々進む。
【００４８】
ステップＳ１０４で異材料間の接触があると判断された場合、材料状態計算手段２１は、以下の数７で示される運動方程式を離散化した数８の式を利用して、粒子Ｉの加速度を求め、これをステップＳ１０３で求めた粒子Ｉの加速度に加算する。材料状態計算手段２１は、粒子Ｉの粒子番号とその加速度とを関連づけて主記憶手段１２に一時保存する（ステップＳ１０５）。ここで、ａｉＩは粒子Ｉの加速度成分、σｉｊＩ,nは時間増分がｎ回進んだ時点の粒子Ｉの応力成分である。この処理が終了したらステップＳ１０６に進む。
【００４９】
【数７】
2014137747.tif 000009
【００５０】
【数８】
2014137747.tif 000010
【００５１】
続いて材料状態計算手段２１は、粒子の速度計算の処理を行う（ステップＳ１０６）。時間増分がｎ回進んだ時点での粒子の速度をｖｉＩ，ｎとすると、時間増分がΔｔ進んだ現時間での加速度をｖｉＩ，ｎ＋１は、ステップＳ１０３および１０５で求められた現時間での粒子Ｉの速度ａｉＩ，ｎ＋１を用いて、以下の数９で示すように求められる。ここで、速度や加速度に限らず、他の変数についてもこのｖやａと同様に時間を表記するものとする。材料状態計算手段２１は、算出されたこの速度についても、粒子Ｉの粒子番号と速度とを関連づけて主記憶手段１２に一時保存する。
【００５２】
【数９】
2014137747.tif 000011
【００５３】
材料状態計算手段２１は、ステップＳ１０６で算出した粒子Ｉの速度から、この粒子と近傍の粒子との間の速度勾配を計算する（ステップＳ１０７）。この速度勾配は、以下の数１０で示すように求められる。
【００５４】
【数１０】
2014137747.tif 000012
【００５５】
そして、材料状態計算手段２１は、ステップＳ１０７で算出した粒子Ｉの速度勾配から、この粒子Ｉの歪速度を計算する（ステップＳ１０８）。粒子Ｉの歪速度は、以下の数１１で示すように求められる。また、これに時間増分Δｔを乗算すれば、数１２で示される現時間の歪増分ΔεｉｊＩ，ｎ＋１を求めることもできる。材料状態計算手段２１による処理はここまでである。
【００５６】
【数１１】
2014137747.tif 000013
【００５７】
【数１２】
2014137747.tif 000014
【００５８】
以上の処理に引き続いて、材料状態判定手段２２は、主記憶手段１２に記憶された材料状態識別情報２６を読み出す（ステップＳ１０９）。材料状態識別情報２６は、前述のステップＳ１０１、および後述のステップＳ１１４とＳ１１７とで主記憶手段１２に記憶されるものであり、当該材料が弾性状態であるか塑性状態であるかを示す二値データを含むものである。
【００５９】
そして材料状態判定手段２２は、読み出した材料状態識別情報２６が弾性状態であるか塑性状態であるかを判断する（ステップＳ１１０）。「弾性状態」であればステップＳ１１１に、「塑性状態」であればステップＳ１１５に、各々処理が進む。
【００６０】
ステップＳ１１０で「弾性状態」と判定された場合、降伏判定手段２３はまず、粒子の弾性の応力計算を行う（ステップＳ１１１）。この計算では、ステップＳ１０８で算出した粒子Ｉの歪増分のベクトル｛ΔεｉｊＩ，ｎ＋１｝を弾性マトリクス［Ｄｅ］に乗算して、以下の数１３に示す粒子Ｉの応力増分のベクトル｛ΔσｉｊＩ，ｎ＋１｝を算出する。さらに、この応力増分ΔσｉｊＩ，ｎ＋１を時間増分Δｔ前の応力σｉｊＩ，ｎに加算し、数１４に示す現時間の応力σｉｊＩ，ｎ＋１を算出する。
【００６１】
【数１３】
2014137747.tif 000015
【００６２】
【数１４】
2014137747.tif 000016
【００６３】
ここで、弾性マトリクス［Ｄｅ］は、以下の数１５にように表される。また、Ｅは粒子Ｉのヤング率、νは粒子Ｉのポアソン比である。
【００６４】
【数１５】
2014137747.tif 000017
【００６５】
そして降伏判定手段２３が粒子の降伏判定を行い（ステップＳ１１２）、これに引き続いて応力修正計算手段２４が粒子の応力修正を行う（ステップＳ１１３）。図５は、図３のステップＳ１１２〜１１３として示した降伏判定手段２３および応力修正計算手段２４による処理の詳細について示すフローチャートである。以下、この図５によって、これらの処理について説明する。
【００６６】
降伏判定手段２３は、まず粒子の静水圧歪増分を算出する（ステップＳ２０１）。この静水圧歪増分ΔεｍＩ，ｎ＋１は、前述のステップＳ１０８で算出された粒子Ｉの歪増分ΔεｉｊＩ，ｎ＋１を用いて、以下の数１６に示されるように算出される。
【００６７】
【数１６】
2014137747.tif 000018
【００６８】
続いて降伏判定手段２３は、前述のステップＳ１０８で算出された現時間の歪増分ΔεｉｊＩ，ｎ＋１から、ステップＳ２０１で算出された現時間の静水圧歪増分ΔεｍＩ，ｎ＋１を減算して、以下の数１７で示される現時間の偏差歪増分ΔｅｉｊＩ，ｎ＋１を算出する（ステップＳ２０２）。
【００６９】
【数１７】
2014137747.tif 000019
【００７０】
続いて降伏判定手段２３は、粒子の時間増分Δｔ前の静水圧応力ΔσｍＩ，ｎを、前述のステップＳ１１１で算出された現時間の応力を用いて、以下の数１８で示されるように算出する（ステップＳ２０３）。
【００７１】
【数１８】
2014137747.tif 000020
続いて降伏判定手段２３は、粒子の時間増分Δｔ前の偏差応力ｓｉｊＩ，ｎを、応力σｉｊＩ，ｎからステップＳ２０３で算出された静水圧応力ΔσｍＩ，ｎを減算して、以下の数１９で示されるように算出する（ステップＳ２０４）。
【００７２】
【数１９】
2014137747.tif 000021
【００７３】
続いて降伏判定手段２３は、粒子の時間増分Δｔ前の偏差応力ｓｉｊＩ，ｎに、ステップＳ２０２で算出された偏差歪み増分に横弾性係数Ｇを乗算して求められる試行偏差応力増分２ＧΔｅｉｊＩ，ｎ＋１を加算して、以下の数２０で示される現時間の試行偏差応力（ｓｉｊＩ，ｎ＋１）Ｔを算出する（ステップＳ２０５）。
【００７４】
【数２０】
2014137747.tif 000022
【００７５】
そして降伏判定手段２３は、ステップＳ２０５で算出された試行偏差応力から、以下の数２１で示される試行相当応力を算出する（ステップＳ２０６）。
【００７６】
【数２１】
2014137747.tif 000023
【００７７】
降伏判定手段２３は、ここまでの計算結果に基づいて、ステップＳ２０６で算出された粒子Ｉの現時間の試行相当応力が降伏応力Ｙ以上であるか否かの判断、即ち粒子Ｉが降伏したか否かの判断を行う（ステップＳ２０７）。この試行相当応力が降伏応力Ｙ未満であれば、粒子Ｉは降伏していないので「弾性状態」であると判断できる。この場合は以下のステップＳ２０８〜２１３に示す応力修正（図３のステップＳ１１３）を省略して図３のステップＳ１１４へ進む。
【００７８】
試行相当応力が降伏応力Ｙ以上となった場合（ステップＳ２０７がイエス）、この場合は粒子Ｉが降伏して「弾性状態」となったと判断されるので、応力修正計算手段２４が以下の数２２で示される粒子Ｉの相当塑性歪増分を計算する（ステップＳ２０８）。ここでＨは、相当塑性歪と相当応力との間の関係を示す以下の数２３の関数である。
【００７９】
【数２２】
2014137747.tif 000024
【００８０】
【数２３】
2014137747.tif 000025
【００８１】
ここで、数２３の関数Ｈが非線形な関数である場合には、この相当塑性歪増分を求めるには非線形方程式の反復法、例えばニュートン法を用いることができる。Ｈが線形な関数である場合には、この相当塑性歪増分は容易に求めることができる。
【００８２】
そして、応力修正計算手段２４は、ステップＳ２０８で求められた相当塑性歪増分を利用して、以下の数２４で示される粒子Ｉの相当応力σを計算する（ステップＳ２０９）。ここで、モデルデータ２７で仮定された材料の塑性が等方性硬化であると登録されている場合には、これを降伏した粒子Ｉの相当応力を時間増分Δｔ後に用いる降伏応力として、粒子Ｉの粒子番号と関連付けて主記憶手段１２に一時的に保存する。
【００８３】
【数２４】
2014137747.tif 000026
【００８４】
これに引き続いて応力修正計算手段２４は、ステップＳ２０９で求められた相当応力σを利用して、以下の数２５で示される粒子Ｉの現時間の偏差応力ｓｉｊＩ，ｎ＋１を計算する（ステップＳ２１０）。
【００８５】
【数２５】
2014137747.tif 000027
【００８６】
さらに応力修正計算手段２４は、ステップＳ２０１で求められた静水圧歪増分ΔεｍＩ，ｎ＋１を利用して、以下の数２６で示される粒子Ｉの現時間の静水圧応力増分ΔσｍＩ，ｎ＋１を算出する（ステップＳ２１１）。ここで、Ｋは体積弾性率である。
【００８７】
【数２６】
2014137747.tif 000028
【００８８】
引き続いて応力修正計算手段２４は、ステップＳ２１１で求められた静水圧応力増分ΔσｍＩ，ｎ＋１を利用して、以下の数２７で示される粒子Ｉの現時間の静水圧応力σｍＩ，ｎ＋１を算出する（ステップＳ２１２）。
【００８９】
【数２７】
2014137747.tif 000029
そして応力修正計算手段２４は、粒子Ｉの現時間の偏差応力に静水圧応力を加算して、以下の数２８で示される粒子Ｉの応力を算出する（ステップＳ２１２）。
【００９０】
【数２８】
2014137747.tif 000030
【００９１】
図３に戻って、応力修正計算手段２４は、材料状態識別情報２６を更新する（ステップＳ１１４）。Ｓ１１２〜１１３（図５）で示した処理の結果に従い、粒子Ｉが降伏していれば「塑性状態」、降伏していなければ「弾性状態」を粒子番号と対応づけて保存する。
【００９２】
一方、前述のステップＳ１１０で「塑性状態」と判定された場合、降伏判定手段２３は、粒子Ｉの弾塑性の応力を計算する（ステップＳ１１５）。より具体的には、以下の数２９によって、ステップＳ１０８で求められた現時間の歪増分のベクトル｛ΔεｉｊＩ，ｎ＋１｝を弾塑性マトリクス［Ｄｅｐ］に乗算して、粒子Ｉの応力増分のベクトル｛ΔσｉｊＩ，ｎ＋１｝を算出する。さらに、ステップＳ１１１（前述の数１４）と同様に、求められた応力増分ΔσｉｊＩ，ｎ＋１を時間増分Δｔ前の応力σｉｊＩ，ｎに加算し、現時間の応力σｉｊＩ，ｎ＋１を算出する。
【００９３】
【数２９】
2014137747.tif 000031
【００９４】
ここで、弾塑性マトリクス［Ｄｅｐ］は以下の数３０で定義される。［Ｄｐ］は対象行列、Ｇは横弾性係数、σは相当応力、Ｈ’は数２３で示される関数Ｈの微分値、ｓｉｉは偏差応力成分、τｉｊは剪断応力成分を表す。
【００９５】
【数３０】
2014137747.tif 000032
【００９６】
そして降伏判定手段２３は、粒子の降伏判定を行う（ステップＳ１１６）。ステップＳ１１５で求められた現時間の現時間の相当応力σが降伏応力Ｙ以上であるか否かの判断、即ち粒子Ｉが降伏したか否かの判断を行う。この相当応力が降伏応力Ｙ未満であれば粒子Ｉは降伏していないので「弾性状態」、降伏応力Ｙ以上であれば粒子Ｉは降伏しているので「塑性状態」であると判断できる。さらに、モデルデータ２７で仮定された材料の塑性が等方性硬化であると登録されている場合には、これを降伏した粒子Ｉの相当応力を時間増分Δｔ後に用いる降伏応力として、粒子Ｉの粒子番号と関連付けて主記憶手段１２に一時的に保存する。
【００９７】
さらに降伏判定手段２３は、前述のステップＳ１１４と同様に、粒子Ｉが降伏していれば「塑性状態」、降伏していなければ「弾性状態」を粒子番号と対応づけて材料状態識別情報２６を更新して保存する（ステップＳ１１７）。
【００９８】
以上ここまでの処理によって更新された材料状態識別情報２６を、終了判定処理手段２５が固定記憶手段１３上のモデルデータ２７に反映させて保存し、また出力手段１５に描画して出力する（ステップＳ１１８）。
【００９９】
そして終了判定処理手段２５は、処理対象となる全ての粒子に対して処理が終了したか否かを判断し（ステップＳ１１９）、終了していなければ次の粒子に進んで（ステップＳ１２０）ステップＳ１０２に戻って処理を繰り返す。全ての粒子に対して処理が終了していれば、処理対象となる全ての時間について終了したか否かを判断し（ステップＳ１２１）、終了していなければ時間増分を加算して次の時間に進んで（ステップＳ１２２）ステップＳ１０２に戻って処理を繰り返す。
【０１００】
処理対象となる全ての粒子および時間に対して処理が終了すれば、処理終了となる。
【０１０１】
（計算の例）
以上ここまでで説明した動作によって、実際の粒子モデルに対して行う計算の例を以下に示す。
【０１０２】
図６は、図２〜４に示した構造解析装置１０で行う計算の対象となる粒子モデル３００について示す説明図である。粒子モデル３００は、直方体３０１および板３０２からなり、粒子の大きさ（直径）を０．２５ｍｍとしてモデル化したものである。
【０１０３】
直方体３０１は、アルミニウムを材料とし、その材料特性はヤング率７１００ｋｇｆ／ｍｍ２（６９．６Ｇｐａ）、ポアソン比０．３３、密度２．７５２Ｅ−１０ｋｇｆ・ｓｅｃ／ｍｍ４（２．７ｇ／ｃｍ３）、塑性構成則（相当塑性歪と相当応力の関係）が以下の数３１のルードウィック則を満たす弾塑性体である。
【０１０４】
【数３１】
2014137747.tif 000033
【０１０５】
板３０２は、ステンレス鋼を材料とし、その材料特性はヤング率２１０００ｋｇｆ／ｍｍ２（２０５．８Ｇｐａ）、ポアソン比０．３３、密度７．９５４Ｅ−１０ｋｇｆ・ｓｅｃ／ｍｍ４（７．８ｇ／ｃｍ３）の剛体である。
【０１０６】
以上の条件で、初速４５００ｍｍ／ｓｅｃ、重力加速度を課した直方体３０１が０．１３５ｍｍ落下して（直方体の１ｍ自由落下を想定）、全方位拘束された板３０２に衝突するという解析を行った。
【０１０７】
図７は、図６に示した粒子モデル３００で、板３０２に衝突する直方体３０１の接触面の中央部に位置する粒子に関して、相当塑性歪と相当応力の関係を図２〜４に示した構造解析装置１０でシミュレートした結果について示すグラフ４００である。このグラフ４００は、塑性構成則の理論値４０１、応力修正を行った本実施形態の粒子法による解析結果４０２、応力修正を行わない従来の粒子法による解析結果４０３、という３本の曲線を含む。
【０１０８】
この図７からわかるように、従来の粒子法による解析結果４０３は、塑性構成則の理論値４０１と比べて明らかに各相当塑性歪に対応する相当応力が小さくなり、理論値と比べての最大誤差は約４２％程度となる。
【０１０９】
これに対して本実施形態の粒子法による解析結果４０２では、塑性構成則の理論値と各相当塑性歪に対応する相当応力は比較的よく一致しており、理論値と比べての最大誤差は約０．９％程度である。従って、本実施形態によれば既存の粒子法と比べて誤差を約４０％低減させることが可能であることを確認できた。
【０１１０】
（実施形態の全体的な動作）
次に、上記の実施形態の全体的な動作について説明する。
本実施形態に係る構造解析方法は、構造物を複数の粒子でモデル化し、各粒子の状態を含むモデルデータを予め記憶し、このモデルデータに対して粒子法による構造解析を行う構造解析装置１０にあって、モデルデータの各粒子に対して材料状態計算手段が加速度を計算し（図３・ステップＳ１０３）、この加速度から時間増分経過後の歪増分を材料状態計算手段が計算し（図３・ステップＳ１０４〜１０８）、記憶された各粒子の状態に材料状態判定手段が歪増分を加算し、時間増分経過後の当該粒子が弾性状態であるか塑性状態であるかを材料状態判定手段が判断し（図３・ステップＳ１１０）、当該粒子の状態が弾性状態である場合に降伏判定手段がその粒子の弾性応力を計算し（図３・ステップＳ１１１）、この弾性応力からその粒子が降伏したか否かを降伏判定手段が判断し（図３・ステップＳ１１２）、粒子が降伏した場合に応力修正計算手段が当該粒子の応力を修正計算して当該粒子の状態を修正記憶し（図３・ステップＳ１１３）、記憶された各粒子の状態を終了判定処理手段が外部に出力する（図３・ステップＳ１１８）。
【０１１１】
ここで、上記各動作ステップについては、これをコンピュータで実行可能にプログラム化し、これらを前記各ステップを直接実行する構造解析装置１０が備えるプロセッサ１１に実行させるようにしてもよい。本プログラムは、非一時的な記録媒体、例えば、ＤＶＤ、ＣＤ、フラッシュメモリ等に記録されてもよい。その場合、本プログラムは、記録媒体からコンピュータによって読み出され、実行される。
この動作により、本実施形態は以下のような効果を奏する。
【０１１２】
本実施形態によれば、構造物の弾塑性特性を考慮する必要がある構造解析に対しても、その構造物の弾性状態から塑性状態への遷移に伴う歪や応力の修正を行いつつ、粒子法を適用した構造解析を行うことができる。これによって、高精度のシミュレーションを短時間で行うことが可能となる。その実際の成果は、図６〜７で示した通りである。
【０１１３】
これまで本発明について図面に示した特定の実施形態をもって説明してきたが、本発明は図面に示した実施形態に限定されるものではなく、本発明の効果を奏する限り、これまで知られたいかなる構成であっても採用することができる。
【０１１４】
上述した実施形態について、その新規な技術内容の要点をまとめると、以下のようになる。なお、上記実施形態の一部または全部は、新規な技術として以下のようにまとめられるが、本発明は必ずしもこれに限定されるものではない。
【０１１５】
（付記１）  構造物を複数の粒子でモデル化し、前記各粒子の状態を含むモデルデータを予め記憶し、このモデルデータに対して粒子法による構造解析を行う構造解析装置であって、
前記モデルデータの前記各粒子に対して加速度を計算し、この加速度から時間増分経過後の歪増分を計算する材料状態計算手段と、
記憶された前記各粒子の状態と前記歪増分とから、前記時間増分経過後の当該粒子が弾性状態であるか塑性状態であるかを判断する材料状態判定手段と、
当該粒子の状態が弾性状態である場合にその粒子の弾性応力を計算して降伏したか否かを判断する降伏判定手段と、
前記粒子が降伏した場合に当該粒子の応力を修正計算して当該粒子の状態を修正記憶する応力修正計算手段と、
記憶された前記各粒子の状態を外部に出力する終了判定処理手段と
を有することを特徴とする構造解析装置。
【０１１６】
（付記２）  前記終了判定処理手段が、
対象となる全ての粒子および時間に対して処理が終了したか否かを判断し、この処理が終了したと判断された場合に記憶された前記各粒子の状態を外部に出力することを特徴とする、付記１に記載の構造解析装置。
【０１１７】
（付記３）  前記降伏判定手段が、当該粒子の試行相当応力を計算することによってその粒子が降伏したか否かを判断する機能を有すると共に、
前記応力修正計算手段が、前記粒子が降伏した場合に前記試行相当応力に基づいて当該粒子の応力を修正計算する機能を有することを特徴とする、付記１に記載の構造解析装置。
【０１１８】
（付記４）  前記材料状態計算手段が、前記各粒子に対して近傍粒子を探索して異材料間の接触の有無を判定し、接触がある場合に離散化された運動方程式から算出された加速度を当該粒子の加速度に加算する機能を有することを特徴とする、付記１に記載の構造解析装置。
【０１１９】
（付記５）  前記降伏判定手段が、当該粒子の状態が塑性状態である場合に、その粒子の弾塑性応力を計算して降伏したか否かを判断する機能を有することを特徴とする、付記１に記載の構造解析装置。
【０１２０】
（付記６）  構造物を複数の粒子でモデル化し、前記各粒子の状態を含むモデルデータを予め記憶し、このモデルデータに対して粒子法による構造解析を行う構造解析装置にあって、
前記モデルデータの前記各粒子に対して材料状態計算手段が加速度を計算し、
この加速度から時間増分経過後の歪増分を前記材料状態計算手段が計算し、
記憶された前記各粒子の状態に材料状態判定手段が前記歪増分を加算し、
前記時間増分経過後の当該粒子が弾性状態であるか塑性状態であるかを前記材料状態判定手段が判断し、
当該粒子の状態が弾性状態である場合に降伏判定手段がその粒子の弾性応力を計算し、
この弾性応力からその粒子が降伏したか否かを前記降伏判定手段が判断し、
前記粒子が降伏した場合に応力修正計算手段が当該粒子の応力を修正計算して当該粒子の状態を修正記憶し、
記憶された前記各粒子の状態を終了判定処理手段が外部に出力する
ことを特徴とする構造解析方法。
【０１２１】
（付記７）  前記終了判定処理手段前記各粒子の状態を外部に出力する処理が、
対象となる全ての粒子および時間に対して処理が終了したか否かを判断し、この処理が終了したと判断された場合に記憶された前記各粒子の状態を外部に出力することを特徴とする、付記６に記載の構造解析方法。
【０１２２】
（付記８）  前記粒子が降伏したか否かを判断する処理が、当該粒子の試行相当応力を計算することによってその粒子が降伏したか否かを判断するものであると共に、
前記粒子の応力を修正計算する処理が、前記粒子が降伏した場合に前記試行相当応力に基づいて当該粒子の応力を修正計算するものであることを特徴とする、付記６に記載の構造解析方法。
【０１２３】
（付記９）  前記各粒子に対して加速度を計算する処理が、前記各粒子に対して近傍粒子を探索して異材料間の接触の有無を判定し、接触がある場合に離散化された運動方程式から算出された加速度を当該粒子の加速度に加算する処理を含むことを特徴とする、付記６に記載の構造解析方法。
【０１２４】
（付記１０）  前記粒子が降伏したか否かを判断する処理が、当該粒子の状態が塑性状態である場合に、その粒子の弾塑性応力を計算して降伏したか否かを判断するものであることを特徴とする、付記６に記載の構造解析方法。
【０１２５】
（付記１１）  構造物を複数の粒子でモデル化し、前記各粒子の状態を含むモデルデータを予め記憶し、このモデルデータに対して粒子法による構造解析を行う構造解析装置にあって、
前記構造解析装置が備えるプロセッサに、
前記モデルデータの前記各粒子に対して加速度を計算する手順、
この加速度から時間増分経過後の歪増分を計算する手順、
記憶された前記各粒子の状態に前記歪増分を加算する手順、
前記時間増分経過後の当該粒子が弾性状態であるか塑性状態であるかを判断する手順、
当該粒子の状態が弾性状態である場合にその粒子の弾性応力を計算する手順、
この弾性応力からその粒子が降伏したか否かを判断する手順、
前記粒子が降伏した場合に当該粒子の応力を修正計算して当該粒子の状態を修正記憶する手順、
および記憶された前記各粒子の状態を終了判定処理手段が外部に出力する手順
を実行させることを特徴とする構造解析プログラム。
【産業上の利用可能性】
【０１２６】
本発明は構造物の変形を解析する用途に対して幅広く適用可能であるが、特に金属の切削加工やプレス加工のシミュレーションなどのような、構造物の弾性状態から塑性状態への遷移を伴う変形の解析に適している。
【符号の説明】
【０１２７】
１、１０  構造解析装置
２、２７  モデルデータ
３、２１  材料状態計算手段
４、２２  材料状態判定手段
５、２３  降伏判定手段
６、２４  応力修正計算手段
７、２５  終了判定処理手段
１１  プロセッサ
１２  主記憶手段
１３  固定記憶手段
１４  入力手段
１５  出力手段
２６  材料状態識別情報
３００  粒子モデル
３０１  直方体
３０２  板
４００  グラフ
４０１  理論値
４０２、４０３  解析結果
【図１】
2014137747.tif 000034
【図２】
2014137747.tif 000035
【図３】
2014137747.tif 000036
【図４】
2014137747.tif 000037
【図５】
2014137747.tif 000038
【図６】
2014137747.tif 000039
【図７】
2014137747.tif 000040