JP2009211681A - 超弾性材料の構成式の係数算出装置、係数算出方法および係数算出用プログラム - Google Patents

Abstract

【課題】Mooney-Rivlinモデルの式の係数について、応力シミュレーションをより安定して行うために必要な安定条件を満たす係数を算出する装置および方法を提供する。
【解決手段】超弾性材料の構成式に使用されるひずみエネルギー密度関数の一つであるMooney-Rivlinモデルに含まれる係数の値を、材料試験の測定値と係数算出式とから推定する係数算出装置において、測定値および計算過程のデータを記憶する記憶装置３と記憶装置３に接続された計算装置２とで構成され、計算装置２は、制約条件付き非線形最小２乗法により係数の近似解を算出する係数算出手段２２と、近似解の誤差を評価し、最適解の候補を更新する誤差評価・解更新手段２３と、近似解の繰り返し計算において、収束状態を評価する収束評価手段２４と、収束困難な場合、最適解の候補を係数の初期値とし、制約条件付き非線形最小２乗法のパラメータを初期化して、制約条件付き非線形最小２乗法をやり直す初期変更手段２５と、を備える。
【選択図】図１

Description

本発明は超弾性材料のひずみエネルギー密度関数の一つであるMooney-Rivlinモデルの係数を算出する方法について、応力シミュレーション時に安定して解を得られる条件を満たす係数を算出可能な係数算出装置、係数算出方法および係数算出プログラムに関する。

ゴムやエラストマーのような樹脂材料は大きく変形することが可能であり、超弾性材料あるいは超弾性体と呼ばれる。超弾性材料の力学構成式では、ひずみエネルギー密度関数を主ストレッチで偏微分することにより、主ストレッチ方向の応力が得られる。ストレッチとは変形後の長さを元の長さで除した値であり、主応力方向のストレッチが主ストレッチである。主応力の方向は直交する３つの方向があることから、主ストレッチも３つ存在する。また、ひずみエネルギー密度関数はそれら３つの主ストレッチの関数として表される。

前記ひずみエネルギー密度関数には、Mooney-Rivlinモデル、Aruda-Boyceモデルなど複数のモデルが提案されている。ＡＮＳＹＳなどの汎用有限要素解析プログラムでは、それらのモデルのひずみエネルギー密度関数の係数を入力することにより、超弾性材料の応力シミュレーションなど各種のシミュレーションを行うことが可能である。

ひずみエネルギー密度関数の中で、Mooney-Rivlinモデルは実験結果へのフィッティング精度が高く、広く使用されているモデルあり、以下の式（１）（ひずみエネルギー密度関数式）で表される。

以下では、従来の係数算出方法について述べる。

非特許文献１によれば、係数Ａ_ｉ，ｊの算出には、図１３〜１５に示す単軸引張試験、等方２軸引張試験、純せん断試験のいずれか一つ以上の材料試験結果が必要である。これらの３つの試験では、負荷を与える方向およびそれに垂直な方向が主応力方向となり、その方向のストレッチが主ストレッチとなる。材料試験結果としては、変位を与える方向のストレッチλと同方向の応力σが出力される。その応力σは、超弾性材料の力学構成式から下の式（５）で与えられる。

式（５）のＺ_ｉ，ｊ（λ）には係数Ａ_ｉ，ｊが含まれておらず、λのみの関数である。材料試験結果のストレッチλ_ｋとそれに対応する応力をσ_ｋ、前記ストレッチλ_ｋを式（５）に代入して得られる応力をσ（λ_ｋ）とすると、データ数ｎ個の材料試験結果と式（５）との２乗誤差Ｅｒｒは次の式（６）で表される。

式（６）で、ｎは実測のデータ数である。２乗誤差Ｅｒｒが最小となるには、式（６）を係数Ａ_ｉ，ｊで偏微分したものが０になることが必要である。

式（７）より下の式（８）が成り立つ。

また、式（５）から、次の式（９）が成り立つ。

式（７）に式（５）、式（９）を代入すると、次の式（１０）が成り立つ。

全ての係数に対して式（１０）が成り立つことから、式（１０）のｉ＋ｊ、ｐ＋ｑ、を１からＮまでは振ると、以下の方程式（１１）を得る。

式（１１）を解くことにより係数Ａ_ｉ，ｊを得る。また、式（１１）を解くことにより、解を得る解法は線形最小２乗法と呼ばれている。

従来方法では、図１１に示すように、材料試験結果を読み込むための入力装置１０１と、係数Ａ_ｉ，ｊを算出する計算装置１０２と、計算過程の試験データを記憶する記憶装置１０３と、求めた係数を出力する出力装置１０４から構成されている。計算装置１０２は、材料試験結果と次数Ｎを読み込み記憶装置１０３に記憶する算出手段１２１と、線形最小２乗法により係数Ａ_ｉ，ｊを求める係数算出手段１２２とから成る。

このような構成を有する従来の係数算出装置は次のように動作する。

即ち、図１２に示すように、算出手段１２１では、入力装置１０１から入力された試験データのストレッチλ_ｋ、σ_ｋ（ｋ＝１・・・ｎ）、およびMooney-Rivlin式の次数Ｎを記録装置１０３に記録し、線形最小２乗法による係数算出手段１２２により記録装置に入力されたデータから方程式（１１）を作成し、前記方程式を解くことにより係数Ａ_ｉ，ｊを算出する。

一方、非特許文献２によれば、応力シミュレーションの計算を安定に行うために、係数Ａ_ｉ，ｊに以下の制約がある。

従来手法では算出された係数Ａ_ｉ，ｊは、式（１２）の制約を必ずしも満たすとは限らない。また、以降、式（１２）を安定条件と呼ぶ。
SIMULIA社「Abaqus Version６，７ Theory Manual ４．６Large-strain elasticity」 サイバネットシステム株式会社「ANSYS構造非線形セミナー」pp.５−５

しかしながら、上述したような従来手法では以下のような問題点がある。

第１の問題点は、従来手法で算出した係数が必ずしも安定条件を満たさないため、応力シミュレーション中に本来生じない変形が発生することにより計算が収束せず、計算結果を得られないということである。その理由は、用いる試験データのばらつきや、必要な試験方法の欠如による。

第２の問題点は、前記第１の問題を回避するためには理想的に３種類の試験（単純引張試験、等方２軸引張試験、純せん断試験）結果が必要であり、コストが掛かるということである。その理由は、３つの内一つ以上の材料試験結果が欠けると、前記欠けた試験の変形状態が係数に考慮されない。そのため、安定条件を満たすためには理想的に３つの材料試験結果が必要だからである。

本発明の目的は、超弾性材料のひずみエネルギー密度関数の一つであるMooney-Rivlinモデルの係数について、安定条件を満たす係数を算出する方法を提供することにある。

上記の課題を解決するために、請求項１に記載の発明は、超弾性材料の構成式に使用されるひずみエネルギー密度関数の一つであるMooney-Rivlinモデルに含まれる係数の値を、材料試験の測定値と係数算出式とから推定する係数算出装置において、前記測定値および計算過程のデータを記憶する記憶装置と該記憶装置に接続された計算装置とで構成され、前記計算装置は、制約条件付き非線形最小２乗法により係数の近似解を算出する係数算出手段と、前記近似解の誤差を評価し、最適解の候補を更新する誤差評価・解更新手段と、前記近似解の繰り返し計算において、前記近似解の収束状態を評価する収束評価手段と、前記収束評価手段の評価に基づき収束困難な場合、最適解の候補を係数の初期値とし、制約条件付き非線形最小２乗法のパラメータを初期化して、制約条件付き非線形最小２乗法をやり直す初期変更手段と、を備えることを特徴とする。

また、請求項２に記載の発明は、材料試験の測定値および計算過程のデータを記憶する記憶装置と該記憶装置に接続された計算装置とで構成され、Mooney-Rivlinモデルのひずみエネルギー密度関数式に含まれる複数の係数の値を材料試験の測定値から推定する係数算出装置における係数算出方法において、制約条件付き非線形最小２乗法により係数の近似解を算出し、前記近似解の誤差を評価し、最適解の候補を更新し、前記近似解の繰り返し計算において、前記近似解の収束状態を評価し、前記評価に基づき収束困難な場合には最適解の候補を係数の初期値として制約条件付き非線形最小２乗法をやり直すことを特徴とする。

また、請求項３に記載の発明は、材料試験の測定値および計算過程のデータを記憶する記憶装置と該記憶装置に接続された計算装置とで構成され、Mooney-Rivlinモデルのひずみエネルギー密度関数式に含まれる複数の係数の値を材料試験の測定値から推定するコンピュータに、制約条件付き非線形最小２乗法により係数の近似解を算出する係数算出ステップと、前記近似解の誤差を評価し、最適解の候補を更新する誤差評価・解更新ステップと、前記近似解の繰り返し計算において、前記近似解の収束状態を評価する収束評価ステップと、前記収束評価ステップの評価に基づき収束困難な場合、最適解の候補を係数の初期値とし、制約条件付き非線形最小２乗法のパラメータを初期化して、制約条件付き非線形最小２乗法をやり直す初期変更ステップと、を実行させるための係数算出プログラムであることを特徴とする。

また、請求項４に記載の発明は、請求項１に記載の係数算出装置において、前記収束評価手段が、さらに前記制約条件付き非線形最小２乗法を終了させる条件を有し、前記終了条件を満たす場合に、前記制約条件付き非線形最小２乗法を終了させる手段を含むことを特徴とする。

また、請求項５に記載の発明は、請求項２に記載の係数算出方法において、前記収束困難な場合で、且つ前記制約条件付き非線形最小２乗法を終了させる条件を満たす場合に、前記制約条件付き非線形最小２乗法を終了させることを特徴とする。

また、請求項６に記載の発明は、請求項３に記載のプログラムにおいて、前記収束状態を評価するステップが、さらに前記制約条件付き非線形最小２乗法を終了させる条件を有し、前記終了条件を満たす場合に、前記制約条件付き非線形最小２乗法を終了させるステップを実行させるためのプログラムであることを特徴とする。

本発明によれば、制約条件付き非線形最小２乗法により得られる係数の近似解の誤差を評価して最適解の候補を更新し、この近似解の収束状態を評価して収束困難な場合には最適解の候補を係数の初期値として制約条件付き非線形最小２乗法をやり直すことにより、安定条件を満たす解を求めることができるため、材料試験結果にばらつきがあっても、応力シミュレーションをより安定して実施できる。また安定条件を満たす係数を算出するために必要な材料試験を減らすことができ、コスト削減を図ることができる。

以下、本発明の最良の実施形態を添付図面に基づいて説明する。

（第１実施形態）
まず、本発明の第１実施形態にかかる係数算出装置の構成について図１を用いて説明する。

図１を参照すると、本発明の係数算出方法および装置は、キーボードやマウス等の入力装置１と、係数の算出を行う計算装置２と、計算時に計算過程のデータを記憶（保持）するための記憶装置３と、ディスプレイや印刷機あるいはファイル装置等の出力装置４とから構成されている。

計算装置２は、初期値決定手段２１と、制約条件付き非線形最小２乗法による係数算出手段２２と、誤差評価・解更新手段（誤差評価手段）２３と、収束評価手段２４と、初期値変更手段２５と、最適解の決定手段（最適解決定手段）２６とから構成される。また、計算装置２に接続された記憶装置３には、２つの記憶部３１、３２が設けられている。ここで、記憶部３１は、材料試験データやMooney-Rivlinモデルの次数や係数の初期値の記憶部であり、記憶部３２は、最適解の候補の記憶部である。

これらの手段はそれぞれ概略つぎのように動作する。

初期値決定手段２１は、入力装置１によって与えられた数値データを記憶装置３に記憶すると共に、係数Ａ_ｉ，ｊの初期値、制約条件付き非線形最小２乗法で使用するパラメータの初期値を決定する。入力装置１から与えられる数値データは、少なくとも以下の２つである。
１）材料試験結果のストレッチλ_ｋと前記ひずみに対応するσ_ｋ（ｋ＝１・・・ｎ）の組
２）Mooney-Rivlinモデルの次数Ｎ
１）、２）のデータは記憶部３１に記憶される。また２）の次数Ｎから係数Ａ_ｉ，ｊの個数がＮ（Ｎ＋３）／２個に決まり、その初期値を記憶部３１に記憶する。

制約条件付き非線形最小２乗法による係数算出手段２２は、記憶部３１から材料試験データ、次数Ｎ、係数Ａ_ｉ，ｊの初期値を読み出し、係数Ａ_ｉ，ｊの近似解を制約条件付き非線形最小２乗法により算出する。つまり、試験で得られた応力σ_ｋと式（５）から得られる応力の理論値σ_ｋ(λ_ｋ)の誤差が係数Ａ_ｉ，ｊの安定条件内で最小となるよう、係数Ａ_ｉ，ｊの近似解を算出する。制約条件付き非線形最小２乗法の代表的な手法には、ペナルティ法や拡張Lagrange法など各種方法があり、近似解の探索にはGauss -Newton法、準Newton法などがあり、任意の手法が利用可能である。ここで、係数算出式の一例としてペナルティ法の１つである外点法を用いた場合、以下の式となる。

但し、σ_ｋ（ｋ=１…ｎ）は材料試験での応力、λ_ｋは前記応力に対応するひずみ、ｒはペナルティパラメータ、Ｍａｘ関数は２つの引数のうち大きい方の値を返す関数である。

誤差評価・解更新手段２３は、係数算出手段２２で得た近似解の誤差と記憶部３２に記憶された解の候補の誤差を比較し、係数算出手段２２で得た近似解の誤差の方が小さい場合、それを記憶部３２に記憶して解の候補を更新する。

収束評価手段２４は、制約条件付き非線形最小２乗法の収束状態、および近似解の探索の収束状態を評価する。近似解の繰り返し計算として、制約条件付き非線形最小２乗法は収束していないが、近似解の探索を継続可能な場合、係数算出手段２２に戻って近似解の探索を継続する。

初期値変更手段２５は、収束評価手段２４で近似解の探索に失敗したような収束困難な場合に、記憶部３２から最適解の候補である現在の最適解を読み取って係数Ａ_ｉ，ｊの初期値とする。さらに、制約条件付き非線形最小２乗法のパラメータを初期化する。その後、現在の最適解を係数Ａ_ｉ，ｊの初期値として制約条件付き非線形最小２乗法をやり直すよう係数算出手段２２に戻る。

最適解決定手段２６は、収束評価手段２４で制約条件付き非線形最小２乗法が収束したと判断された場合、記憶部３２から係数Ａ_ｉ，ｊの最適解を読み取り、出力装置４に出力する。

以上のように、超弾性材料の構成式に使用されるひずみエネルギー密度関数の一つであるMooney-Rivlinモデルに含まれる係数Ａ_ｉ，ｊの値を、材料試験の測定値と係数算出式とから推定する。

次に、図１のブロック図および図２の流れ図を参照して本実施形態の全体動作について詳細に説明する。

まず、入力装置１から与えられた材料試験結果のストレッチλ_ｋと応力σ_ｋ（ｋ＝１・・・ｎ）、次数Ｎは、図２に示すように初期値決定手段２１を通じて、記憶部３１に記憶される（ステップＡ１）。また、係数Ａ_ｉ，ｊの個数をＮ（Ｎ＋３）／２個としてそれらの初期値が記憶部３１に記憶されると共に、最適解の最初の候補として記憶部３２に記憶される。（ステップＡ２）。

次に、制約条件付き非線形最小２乗法により係数Ａ_ｉ，ｊの近似解を算出する（ステップＡ３）。制約条件付き非線形最小２乗法の代表的な手法にはペナルティ法、拡張Lagrange法など各種手法があり任意の手法が利用可能である。また、近似解の代表的な算出方法には、Newton法、Gauss-Newton法、準Newton法など各種方法があり任意の手法が利用可能である。前記近似解について、安定条件を満たすか（ステップＡ４）、あるいは前記近似解の誤差が記憶部３２に記憶されている最適解の候補より誤差が小さいか（ステップＡ５）を確認する。安定条件を満たさない（ステップＡ４）、あるいは現在の最適解より誤差が大きい（ステップＡ５）場合は、近似解の検索を継続可能か確認する（ステップＡ８）。近似解の探索を継続可能の場合、ステップＡ３へ戻って近似解の探索を継続する。

一方、安定条件を満たし（ステップＡ４）、且つ現在の最適解より誤差が小さい（ステップＡ５）場合は、記憶部３２に保存された最適解を前記近似解に置き換えて更新する（ステップＡ６）。その後、収束状態を確認する（ステップＡ７）。

収束条件を満たさない場合（ステップＡ７）、近似解の探索を継続可能か確認する（ステップＡ８）。近似解の探索を継続可能な場合、ステップＡ３へ戻って近似解の探索を継続する。

ステップＡ８で、近似解の探索を継続できない場合（例えば、オーバーフローを起こしそうになっている、繰り返し計算の上限値を越えているなど）、係数Ａ_ｉ，ｊの値を記憶部３２に記憶された値に変更してそれを初期値し、さらに制約条件付き非線形最小２乗法のパラメータを初期化する（ステップＡ９）。その後、ステップＡ３に戻り、係数Ａ_ｉ，ｊの初期値を現在の最適解として新たに制約条件付き非線形最小２乗法による係数算出を行う。

収束条件を満たす場合（ステップＡ７）、記憶部３２から最適解を読み出して出力装置４に出力する（ステップＡ１０）。

次に、本実施形態の効果について説明する。

本実施形態では、制約条件付き非線形最小２乗法による係数算出手段２２により係数を算出し、誤差評価・解更新手段２３により係数の近似解の誤差を評価して最適解の候補を更新し、収束評価手段２４により近似解の収束状態を評価し、初期値変更手段２５により、収束評価手段２４の評価に基づき収束困難な場合には最適解の候補を係数の初期値として制約条件付き非線形最小２乗法をやり直すことで、超弾性材料のひずみエネルギー密度関数の一つであるMooney-Rivlinモデルの係数について安定条件を満たす解を求めることができる。そのため、材料試験結果にばらつきがあっても、応力シミュレーションをより安定して実施できる。また安定条件を満たす係数を算出するために必要な材料試験を減らすことができ、コスト削減を図ることができる。さらに必要な材料試験に実施不能な試験があるような材料（例えば、純せん断試験が難しいフィルム材）でも安定条件を満たす係数を算出できる。

また、本実施形態では、近似解が求められるたびに誤差を評価して最適解を更新するように構成されているため、収束条件が厳しすぎて収束条件を満足する解が得られなかった場合でも安定条件を満たす次善の解を得ることができる。

（第２実施形態）
次に、本発明の第２実施形態にかかる係数算出装置について図１、図３に基づき説明する。なお、前記第１実施形態と同一または対応する部分には、同一の符号を用いて異なる構成および作用のみを説明する。その他の実施形態および変形例も同様とする。

図３は、本発明の第２実施形態にかかる係数算出装置の動作を示す流れ図である。

図１および図３に示すように、収束評価手段２４が、さらに制約条件付き非線形最小２乗法の終了条件を満たすかを判断する機能を有する。すなわち、収束評価手段２４が、制約条件付き非線形最小２乗法を終了させる条件を有し、終了条件を満たす場合に、制約条件付き非線形最小２乗法を終了させる手段を含んでいる。

次に、図１のブロック図および図３の流れ図を参照して本実施形態の全体動作について詳細に説明するが、ステップＡ８までは、同じ動作であるので、ステップＡ８の後の工程について説明する。

図３に示すように、ステップＡ８で、近似解の探索を継続できない場合（例えば、オーバーフローを起こしそうになっている、繰り返し計算の上限値を越えているなど）、さらに非線形最小２乗法を終了させる条件を評価する（ステップＡ８１）。このステップＡ８１では、前記終了条件（例えば、初期化回数が所定の値を越えた、最後の最適解更新後に複数回の初期化が行われたなど）と比較する。

前記終了条件を満たさない場合、係数Ａ_ｉ，ｊの値を記憶部３２に記憶された値に変更してそれを初期値とし、さらに制約条件付き非線形最小２乗法のパラメータを初期化する（ステップＡ９）。その後、ステップＡ３に戻り、係数Ａ_ｉ，ｊの初期値を現在の最適解として新たに制約条件付き非線形最小２乗法による係数算出を行う。

ステップＡ８１で前記終了条件を満たす場合、ステップＡ１０へ進み、記憶部３２から最適解を読み出して出力装置４に出力する。

このように、本実施形態では、収束評価手段２４で近似解の探索に失敗したような収束困難な場合で、且つ、制約条件付き非線形最小２乗法を終了させる条件を満たす場合に、制約条件付き非線形最小２乗法を終了させる。

次に、本実施形態の効果について説明する。

本実施形態では、第１実施形態の効果に加え、収束困難な場合でもループに陥らずに確実に計算を終了することができ、計算コストの削減を図ることができる。

（第３実施形態）
次に、本発明にかかる第３実施形態について図に基づき説明する。

図４を参照すると、本発明の第３実施形態は、インターネット等のネットワーク５を通じてサーバコンピュータ６とクライアントコンピュータ７とが接続されたネットワークシステムにおいて、クライアントコンピュータ７から係数算出要求をネットワーク５を通じてサーバコンピュータ６へ送信し、この要求を受けたサーバコンピュータ６が係数算出処理を行って求めた係数の値をネットワーク５を通じて要求元のクライアントコンピュータ７へ返却するサービスを有料または無料で提供する。

サーバコンピュータ６には、図１の計算装置２と記憶装置３とが設けられ、クライアントコンピュータ７には、図１の入力装置１と出力装置４とが設けられ、クライアントコンピュータ７に設けられた入力装置１および出力装置４とサーバコンピュータ６に設けられた計算装置２とがネットワーク５を通じて通信可能に接続される。

次に本実施形態の動作を説明する。

クライアントコンピュータ７のユーザは、ネットワーク５を通じてサーバコンピュータ６にクライアントコンピュータ７を接続し、材料試験結果λ_ｋ、σ_ｋ（ｋ＝１・・・ｎ）、次数Ｎをサーバコンピュータ６へ送信する。

サーバコンピュータ６は、係数算出要求を受信すると、それに含まれる材料試験結果λ_ｋ、σ_ｋ（ｋ＝１・・・ｎ）、次数Ｎを記憶装置に記憶し、第１実施形態及び第２実施形態における計算装置２と同様の処理を実行して係数Ａ_ｉ，ｊを算出し、それらの値を含む応答メッセージをネットワーク５を通じてクライアントコンピュータ７へ送信する。

クライアントコンピュータ７は、サーバコンピュータ６から応答メッセージを受信すると、それに含まれる係数Ａ_ｉ，ｊの値を記憶装置に記憶する。記憶された係数Ａ_ｉ，ｊの値は、各種のシミュレーションに利用される。

本実施形態によれば、ＡＳＰ（Application Service Provider）の形式で、超弾性材料のひずみエネルギー密度関数の一つであるMooney-Rivlinモデルの係数について、安定条件を満たす係数を算出するサービスを提供できる。

次に、具体的な実施例を用いて本発明を実施するための最良の形態の動作を説明する。

まず、第１実施形態に対応する実施例１について図２を参照しながら説明する。

図５に、材料試験の対象としてエラストマー樹脂を単軸引張試験した結果を示す。なお、Mooney-Rivlinモデルの次数は３とした。

初期値決定手段２１は、材料試験結果の５４個のストレッチλ_ｋと対応する応力σ_ｋ、および次数３を読み込み、記憶部３１に記憶する（ステップＡ１）。次に３×（３+３）/２=９個の係数Ａ_１，０、Ａ_０，１、Ａ_２，０、Ａ_１，１、Ａ_０，２、Ａ_３，０、Ａ_２，１、Ａ_１，２、Ａ_０，３を生成し、Ａ_０，３を−１に、他の全ての係数に１を代入した後、前記９つの係数の値を最初の最適解として記憶部３２に記憶する（ステップＡ２）。

制約条件付き非線形最小２乗法による係数算出手段２２では、制約条件付き非線形最小２乗法によって係数Ａ_ｉ，ｊの近似解を逐次算出する（ステップＡ３）。本実施例では、ペナルティ法の中でも外点法を使用し、近似解の探索に準Newton法を用いた。前記手法では、近似解の探索回数の上限値を２００、収束条件を式（６）で表される２乗誤差が８０以下とした。

ステップＡ３により、１回目の探索で得た近似解は
（Ａ１，０、Ａ_０，１、Ａ_２，０、Ａ_１，１、Ａ_０，２、Ａ_３，０、Ａ_２，１、Ａ_１，２、Ａ_０，３）
＝（841.989053766045, 707.691421842189, - 18.3673423705181, 59.4131914163343, 95.3784826320926, -620.524096637151, -380.914550813533, -227.762760412208, -133.120143621025）である。前記近似解では、Ａ_３，０＝-620.524096637151であり安定条件を満たしていない（ステップＡ４）。また、探索回数１回は上限の２００回より小さい（ステップＡ８）。そのためステップＡ３へ戻って近似解の検索を継続する。

上記手順を繰り返し、探索回数１２回目で得た近似解は
（Ａ１，０、Ａ_０，１、Ａ_２，０、Ａ_１，１、Ａ_０，２、Ａ_３，０、Ａ_２，１、Ａ_１，２、Ａ_０，３）
＝（-1079.05584357833, 1146.62019191279, - 1269.28541919981, 956.316783202344, 1435.30475104625, 1019.32880953864, -3617.33844981249, 5237.51445012979, -2612.98663236558）である。前記１２回目の近似解は安定条件を満たすため（ステップＡ４）、現在の最適解である係数Ａ_ｉ，ｊの初期値
（Ａ１，０、Ａ_０，１、Ａ_２，０、Ａ_１，１、Ａ_０，２、Ａ_３，０、Ａ_２，１、Ａ_１，２、Ａ_０，３）＝（1,1,1,1,1,1,1,1,-1）の誤差と前記１２回目の近似解の誤差を比較する（ステップＡ５）。現在の最適解の誤差７．１３×１０４よりも、前記１２回目の近似解の誤差５．６１×１０４の方が小さいため、記憶部３２に記録された最適解を前記１２回目の近似解に更新する（ステップＡ６）。前記１２回目の近似解の誤差は、５．６１×１０４であり、収束条件（誤差８０以下）を満たしていない（ステップＡ７）。また、探索回数は１２回であり、探索回数の上限値よりも小さい（ステップＡ８）。そのため、ステップＡ３に戻って近似解の検索を継続する。

上記手順を繰り返したが、計算は一度も収束条件を満たすことなく検索回数の上限２００回に達した。近似解の探索継続は不能と判断される（ステップＡ８）。そのため、記憶部３２から現在の最適解を読み込んで係数Ａ_ｉ，ｊの初期値に設定する。即ち、記憶部３２に保存されている、検索回数１００回目の近似解
（Ａ_１，０、Ａ_０，１、Ａ_２，０、Ａ_１，１、Ａ_０，２、Ａ_３，０、Ａ_２，１、Ａ_１，２、Ａ_０，３）
＝（-1502.08720027842, 1616.10473943169, 14522.3344667581, - 35222.8065812693, 22503.5981350251, 2986.05785736438, -12875.8623480809, 14512.0854090154, -4622.28086978699）を係数Ａ_ｉ，ｊの初期値に設定する。さらに検索回数を０に戻すなど、非線形最小２乗法のパラメータも初期化する（ステップＡ９）。その後、ステップＡ３に戻って制約条件付き非線形最小２乗法をやり直す。

以上の手順を繰り返した結果、ステップＡ９での初期化回数２回、近似解の探索回数７６回で近似解
（Ａ_１，０、Ａ_０，１、Ａ_２，０、Ａ_１，１、Ａ_０，２、Ａ_３，０、Ａ_２，１、Ａ_１，２、Ａ_０，３）
＝（-1461.27569609874, 1573.81842392817, 15313.2966433464, - 36701.50877382, 23133.208695068, 2139.21216686939, -9083.0281259646, 8605.59924755559, -1661.72246506214）を得た。前記近似解は、安定条件を満たしている（ステップＡ４）。また、２乗誤差は８０以下であるため、記憶部３２に記憶された最適解の候補よりも誤差が小さく（ステップＡ５）、記憶部３２の最適解は前記近似解に更新される（ステップＡ６）。さらに収束条件である２乗誤差８０以下を満たす（ステップＡ７）ことから、前記近似解が最終的な解となり、最適解として出力装置４に出力される（ステップＡ１０）。

本発明の手法で得た解と、従来手法で得た解をあわせて図６に示す。また、従来手法および本発明の手法で得た係数を用いた際の応力−ストレッチ曲線と材料試験結果との比較を図７に示す。図７から、本手法で得た係数は実用上十分な精度を有することが分かる。

次に、第２実施形態に対応する実施例２について図３を参照しながら説明する。

図８に、材料試験の対象としてエラストマー樹脂を単軸引張試験した結果を示す。なお、Mooney-Rivlinモデルの次数は実施例１と同様に３とした。

初期値決定手段２１は、試験結果の４２個のストレッチλ_ｋと対応する応力σ_ｋ、および次数３を読み込み、記憶部３１に記憶する（ステップＡ１）。次に３×（３+３）/２=９個の係数Ａ_１，０、Ａ_０，１、Ａ_２，０、Ａ_１，１、Ａ_０，２、Ａ_３，０、Ａ_２，１、Ａ_１，２、Ａ_０，３を生成し、Ａ_０，３を−１に、他の全ての係数に１を代入した後、前記９つの係数の値を最初の最適解として記憶部３２に記憶する（ステップＡ２）。

制約付き非線形最小２乗法による係数算出手段２２では、制約付き非線形最小２乗法によって係数Ａ_ｉ，ｊの近似解を逐次算出する（ステップＡ３）。本実施例では、ペナルティ法の中でも外点法を使用し、近似解の探索に準Newton法を用いた。前記手法では、近似解の探索回数の上限値を２００、収束条件を式（６）で表される２乗誤差が１００以下、制約条件付き非線形最小２乗法の終了条件として初期化回数が４回以上とした。

ステップＡ３により、１回目の探索で得た近似解は
（Ａ_１，０、Ａ_０，１、Ａ_２，０、Ａ_１，１、Ａ_０，２、Ａ_３，０、Ａ_２，１、Ａ_１，２、Ａ_０，３）
＝（-2371.04378286413, -1044.45236127576, -10336.7992923788, -5948.11100572578, -3371.01050198821, -27665.9822035386, -16588.7944256015, -9902.09260413247, -5881.03772830337）である。

前記近似解では、 Ａ_１，０ + Ａ_０，１＜ ０ であり安定条件を満たしていない（ステップＡ４）。また、探索回数１回は上限の２００回より小さい（ステップＡ８）。そのためステップＡ３へ戻って近似解の検索を継続する。

上記手順を繰り返し、探索回数１１回目で得た近似解は
（Ａ_１，０、Ａ_０，１、Ａ_２，０、Ａ_１，１、Ａ_０，２、Ａ_３，０、Ａ_２，１、Ａ_１，２、Ａ_０，３）
＝（-427.968377538411, 523.729618343615, -96.7671977604242, 151.126754460605, 82.3185633774202, 136.195625149819, -763.68533893447, 1356.02136007691, -719.025598788785）である。前記１１回目の近似解は安定条件を満たすため（ステップＡ４）、現在の最適解である係数Ａ_ｉ，ｊの初期値
（Ａ_１，０、Ａ_０，１、Ａ_２，０、Ａ_１，１、Ａ_０，２、Ａ_３，０、Ａ_２，１、Ａ_１，２、Ａ_０，３）＝（1,1,1,1,1,1,1,1,-1）の誤差と前記１１回目の近似解の誤差を比較する（ステップＡ５）。現在の最適解の誤差２．８２×１０^５よりも、前記１１回目の近似解の誤差２．３０×１０^４の方が小さいため、記憶部３２に記録された最適解を前記１１回目の近似解に更新する（ステップＡ６）。前記１１回目の近似解の誤差は、２．３０×１０^４であり、収束条件（誤差１００以下）を満たしていない（ステップＡ７）。また、探索回数は１１回であり、探索回数の上限値よりも小さい（ステップＡ８）。そのため、ステップＡ３に戻って近似解の検索を継続する。

上記手順を繰り返したが、計算は一度も収束条件を満たすことなく検索回数の上限２００回に 達した。近似解の探索継続は不能と判断される（ステップＡ８）。また、初期化回数は０回のため、制約条件付き非線形最小２乗法を継続可能である（ステップ Ａ８１）。そのため、記憶部３２から現在の最適解を読み込んで係数Ａ_ｉ，ｊの初期値に設定し、初期化回数を１増加する。即ち、記憶部３２に保存されている、検索回数１８８回目の近似解
（Ａ_１，０、Ａ_０，１、Ａ_２，０、Ａ_１，１、Ａ_０，２、Ａ_３，０、Ａ_２，１、Ａ_１，２、Ａ_０，３）
＝（-1505.80848438853, 1633.94985908882, 9121.95764249245, -22570.6375043711, 15009.6200366711, 607.051496951451, -2923.64330278058, 2316.5926437936, -0.000156251014545927）を係数Ａ_ｉ，ｊの初期値に設定する。さらに検索回数を０に戻すなど、非線形最小２乗法のパラメータも初期化する（ステップ９）。その後、ステップＡ３に戻って制約付き非線形最小２乗法をやり直す。

以上の手順を繰り返した結果、初期化回数１回、探索回数７５回（繰り返し計算回数２７５ 回）でエラーが発生し、近似解の探索継続は不能と判断される（ステップＡ８）。また、初期化回数は１回であるため、制約条件付き非線形最小２乗法を継続可能である（ステップＡ８１）。そのため、記憶部３２から現在の最適解を読み込んで係数Ａ_ｉ，ｊの初期値に設定し、初期化回数を１増加する。即ち、記憶部３２に保存されている、検索回数１８８回目の近似解
（Ａ_１，０、Ａ_０，１、Ａ_２，０、Ａ_１，１、Ａ_０，２、Ａ_３，０、Ａ_２，１、Ａ_１，２、Ａ_０，３）
＝（-1505.80848438853, 1633.94985908882, 9121.95764249245, -22570.6375043711, 15009.6200366711, 607.051496951451, -2923.64330278058, 2316.5926437936, -0.000156251014545927）を係数Ａ_ｉ，ｊの初期値に設定する。さらに検索回数を０に戻すなど、非線形最小２乗法のパラメータも初期化する（ステップ９）。その後、ステップＡ３に戻って制約付き非線形最小２乗法をやり直す。

その後、繰り返し計算回数４１３回（初期化回数２回）、５５１回（初期化回数３回）、６８９ 回（初期化回数４回）でエラーが発生した。繰り返し計算回数６８９回のエラーでは、初期化回数４回となり制約条件付き非線形最小２乗法の継続が不能と判断 される（ステップＡ８１）。制約条件付き非線形最小２乗法を終了し、最適解の決定を行う（ステップＡ１０）。

ステップＡ１０では、記憶部３２の最適解を読み取り、前記近似解を最終的な最適解として出力装置４に出力される（ステップＡ１０）。出力される最適解は、初期化回数０回、探索回数１８８回の近似解
（Ａ_１，０、Ａ_０，１、Ａ_２，０、Ａ_１，１、Ａ_０，２、Ａ_３，０、Ａ_２，１、Ａ_１，２、Ａ_０，３）
＝（-1505.80848438853, 1633.94985908882, 9121.95764249245, -22570.6375043711, 15009.6200366711, 607.051496951451, -2923.64330278058, 2316.5926437936, -0.000156251014545927）である。

本発明の手法で得た解と、従来手法で得た解をあわせて図９に示す。また、従来手法および本発明の手法で得た係数を用いた際の応力−ストレッチ曲線と材料試験結果との比較を図１０に示す。図１０から、本手法で得た係数は実用上十分な精度を有することが分かる。

以上、本発明の実施の形態および実施例について説明したが、本発明は以上の例に限定されず、その他各種の付加変更が可能である。例えば、係数Ａ_ｉ，ｊの初期値は
（Ａ_１，０、Ａ_０，１、Ａ_２，０、Ａ_１，１、Ａ_０，２、Ａ_３，０、Ａ_２，１、Ａ_１，２、Ａ_０，３）＝（1,1,1,1,1,1,1,1,-1）に限定されるものではない。また、係数Ａ_ｉ，ｊの初期値、制約条件付き非線形最小２乗法で使用するパラメータ、収束条件、検索回数の上限値を入力装置１から入力しても良い。この場合、記憶部３１に前記パラメータが記憶される。 さらに、特に実施例２では、制約条件付き非線形最小２乗法の終了条件として初期化回数４回以上としたが、入力装置１から終了条件となる初期化回数を入力しても良い。この場合、記憶部３１に前記パラメータが記録される。また、制約条件付き非線形最小２乗法の終了条件を、初期化時の解と最適解が一致するか否か（即ち、初期化後に最適解が更新されたか否か）としても良い。ペナルティパラメータにしきい値を設け、それを制約条件付き非線形最小２乗法の終了条件としても良い。繰り返し計算回数に上限を設け、それを制約条件付き非線形最小２乗法の 終了条件としても良い。

また、本発明の実施形態および実施例では完全非圧縮性のMooney-Rivlinモデルの係数算出に本発明を適用したが、微圧縮を仮定したMooney-RivlinモデルやPolynomial Formモデルの体積項以外の係数を求める際にも適用できる。これは、前記モデルの体積項が体積試験（直交する３方向に同一変位を加える試験）のみから算出されるためである。

さらに、本発明は、材料構成式に含まれる、制約がある複数の係数の値を材料試験の測定値から推定する場合に広く適用することが可能である。

本発明によれば、力学構成式の係数に安定条件がある力学モデルの係数を、材料試験の測定値から推定する分野に利用できる。特に、超弾性材料のひずみエネルギー密度関数の一つであるMooney-Rivlinモデルの係数を算出する用途に適用できる。他には、Polynomial Formモデルの係数を算出する用途に適用できる。

本発明の第１実施形態にかかる係数算出装置の構成を示すブロック図である。 図１に示す係数算出装置の動作を示す流れ図である。 本発明の第２実施形態にかかる係数算出装置の動作を示す流れ図である。 本発明の第３実施形態にかかる係数算出装置の構成を示すブロック図である。 実施例１で使用する単軸引張試験結果を示す説明図である。 図５の材料試験結果に対して、第１実施形態の係数算出装置により算出した係数および従来手法で算出した係数を示す説明図である。 図５の材料試験結果と、図６の係数を用いた際の応力−ストレッチ曲線を示すグラフである。 実施例２で使用する単軸引張試験結果を示す説明図である。 図８の材料試験結果に対して、第１実施形態の係数算出装置により算出した係数および従来手法で算出した係数を示す説明図である。 図８の材料試験結果と、図９の係数を用いた際の応力−ストレッチ曲線を示すグラフである。 従来手法の構成を示すブロック図である。 従来手法の動作を示す流れ図である。 本発明の手法で使用可能な材料試験方法の一つで、単軸引張試験の模式図である。 本発明の手法で使用可能な材料試験方法の一つで、等方２軸引張試験の模式図である。 本発明の手法で使用可能な材料試験方法の一つで、純せん断試験の模式図である。

符号の説明

１…入力装置
２…計算装置
２１…初期値決定手段
２２…係数算出手段
２３…誤差評価・解更新手段
２４…収束評価手段
２５…初期値変更手段
２６…最適解決定手段
３…記憶装置
３１…記憶部
３２…記憶部
４…出力装置
５…ネットワーク
６…サーバコンピュータ
７…クライアントコンピュータ

Claims (6)

1. 超弾性材料の構成式に使用されるひずみエネルギー密度関数の一つであるMooney-Rivlinモデルに含まれる係数の値を、材料試験の測定値と係数算出式とから推定する係数算出装置において、
   前記測定値および計算過程のデータを記憶する記憶装置と該記憶装置に接続された計算装置とで構成され、
   前記計算装置は、
   制約条件付き非線形最小２乗法により係数の近似解を算出する係数算出手段と、
   前記近似解の誤差を評価し、最適解の候補を更新する誤差評価・解更新手段と、
   前記近似解の繰り返し計算において、前記近似解の収束状態を評価する収束評価手段と、
   前記収束評価手段の評価に基づき収束困難な場合、最適解の候補を係数の初期値とし、制約条件付き非線形最小２乗法のパラメータを初期化して、制約条件付き非線形最小２乗法をやり直す初期変更手段と、を備えたことを特徴とする係数算出装置。
2. 材料試験の測定値および計算過程のデータを記憶する記憶装置と該記憶装置に接続された計算装置とで構成され、Mooney-Rivlinモデルのひずみエネルギー密度関数式に含まれ
   る複数の係数の値を材料試験の測定値から推定する係数算出装置における係数算出方法において、
   制約条件付き非線形最小２乗法により係数の近似解を算出し、
   前記近似解の誤差を評価し、最適解の候補を更新し、
   前記近似解の繰り返し計算において、前記近似解の収束状態を評価し、
   前記評価に基づき収束困難な場合には最適解の候補を係数の初期値として制約条件付き非線形最小２乗法をやり直すことを特徴とする係数算出方法。
3. 材料試験の測定値および計算過程のデータを記憶する記憶装置と該記憶装置に接続された計算装置とで構成され、Mooney-Rivlinモデルのひずみエネルギー密度関数式に含まれ
   る複数の係数の値を材料試験の測定値から推定するコンピュータに、
   制約条件付き非線形最小２乗法により係数の近似解を算出する係数算出ステップと、
   前記近似解の誤差を評価し、最適解の候補を更新する誤差評価・解更新ステップと、
   前記近似解の繰り返し計算において、前記近似解の収束状態を評価する収束評価ステップと、
   前記収束評価ステップの評価に基づき収束困難な場合、最適解の候補を係数の初期値とし、制約条件付き非線形最小２乗法のパラメータを初期化して、制約条件付き非線形最小２乗法をやり直す初期変更ステップと、を実行させるための係数算出プログラム。
4. 請求項１に記載の係数算出装置において、
   前記収束評価手段が、さらに前記制約条件付き非線形最小２乗法を終了させる条件を有し、前記終了条件を満たす場合に、前記制約条件付き非線形最小２乗法を終了させる手段を含むことを特徴とする係数算出装置。
5. 請求項２に記載の係数算出方法において、
   前記収束困難な場合で、且つ前記制約条件付き非線形最小２乗法を終了させる条件を満たす場合に、前記制約条件付き非線形最小２乗法を終了させる
   ことを特徴とする係数算出方法。
6. 請求項３に記載のプログラムにおいて、
   前記収束状態を評価するステップが、さらに前記制約条件付き非線形最小２乗法を終了させる条件を有し、前記終了条件を満たす場合に、前記制約条件付き非線形最小２乗法を終了させるステップを実行させるためのプログラム。

Images