JP5782604B2

JP5782604B2 - 情報処理装置及びプログラム

Info

Publication number: JP5782604B2
Application number: JP2012197851A
Authority: JP
Inventors: 田中　真人; 真人田中; 竜二表; 崇笹川; 正毅藤川
Original assignee: University of the Ryukyus; Toyota Central R&D Labs Inc
Current assignee: University of the Ryukyus; Toyota Central R&D Labs Inc
Priority date: 2012-09-07
Filing date: 2012-09-07
Publication date: 2015-09-24
Anticipated expiration: 2032-09-07
Also published as: JP2014052909A; WO2014038729A2; EP2893464A2; WO2014038729A3; CN104603771B; CN104603771A; US20150205896A1

Description

本発明は、情報処理装置及びプログラムに係り、特に、関数の微分計算を行う情報処理装置及びプログラムに関する。

近年、高機能な汎用有限要素法（ＦｉｎｉｔｅＥｌｅｍｅｎｔＭｅｔｈｏｄ、以下ＦＥＭと略記）解析ソフトウェアが数多く市販されてきており、これらの汎用ソフトウェアを利用して効率よく設計業務を進めることは、ものづくりの現場において、ごく一般的に行われている。しかしながら、ユーザが抱える解析業務の中には、これらの汎用ソフトウェアが備える機能の範囲を超えた特殊な解析技術を要求されることがしばしばある。

このような問題に対処するため、多くの汎用ＦＥＭ解析ソフトウェアはユーザサブルーチン機能を提供しており、ユーザ自身がカスタマイズして独自の解析技術やモデルを汎用ソフトウェアの中に組み込むことを可能にしている。通常、汎用ＦＥＭソフトウェアにおける材料構成則のユーザサブルーチンでは、所望の材料構成則を実装するために、与えられた変位・ひずみ量に対して、応力値と接線剛性を求める際に必要となる応力ひずみマトリクス（材料ヤコビアンと呼ぶ）を計算してメインプログラムへ返すことが要求される。接線剛性及び材料ヤコビアンは、Ｎｅｗｔｏｎ−Ｒａｐｈｓｏn型の収束計算で必要であり、本来は応力増分アルゴリズムと整合した値を返さねばならない。

特に、時間増分を大きく取りたい場合や、材料非線形問題・大変形問題のような非線形性の強い問題へ適用する場合などでは、正確な整合接線剛性及び材料ヤコビアンの計算値が不可欠である。さらに、整合接線剛性はＮｅｗｔｏｎ−Ｒａｐｈｓｏｎ法の２次収束性のみならず、正確な感度や座屈固有値を得るためにも重要である。しかし、材料構成則が複雑であるほどに解析的に求めることは難しく、計算を一部でも間違えると、最悪の場合解が発散してしまうこともあるため、その計算には細心の注意が要求される。また、材料構成則によっては現実的に導出自体が不可能なものも少なくない。

材料ヤコビアンは、応力をひずみで微分することで得られる。非特許文献１、２では、この材料ヤコビアンの複雑な解析解の導出を省略するため、以下の（１）式の前方オイラー法を用いた数値微分を利用している。

ここで、ｆ（ｘ）はスカラー関数、ｆ’（ｘ）は関数ｆ（ｘ）の一階微分、Δｘは微小な変動量である。

一方、非特許文献３では、丸め誤差のない数値微分法として、以下の（２）式の複素数階微分近似（Ｃｏｍｐｌｅｘ−ＳｔｅｐＤｅｒｉｖａｔｉｖｅＡｐｐｒｏｘｉｍａｔｉｏｎ）を提唱しており、その優れた効果を報告している。

ここで、ｉは虚数単位、Ｉｍは虚部を取り出す演算子である。複素数階微分法は、微分操作を複素数平面に拡張することによって、一階微分近似に関してはいかなる小さな変動量Δｘを与えても丸め誤差を決して生じさせないという強力な性質を持つ。複素数階微分法を用いれば、問題に非依存な変動量Δｘの設定が可能となり、汎用的で高精度な微分近似が得られる。非特許文献４では、この複素数階微分法を用いて、上記の非特許文献１、２の手法を拡張し、整合接線剛性の高精度な数値近似手法の導出を行っている。

Ｍｉｃｈｅ，Ｃ．，ＮｕｍｅｒｉｃａｌＣｏｍｐｕｔａｔｉｏｎｏｆａｌｇｏｒｉｔｈｍｉｃ（ｃｏｎｓｉｓｔｅｎｔ）ｔａｎｇｅｎｔｍｏｄｕｌｉｉｎｌａｒｇｅ−ｓｔｒａｉｎｃｏｍｐｕｔａｔｉｏｎａｌｉｎｅｌａｓｔｉｃｉｔｙ，ＣｏｍｐｕｔｅｒＭｅｔｈｏｄｓｉｎＡｐｐｌｉｅｄＭｅｃｈａｎｉｃｓａｎｄＥｎｇｉｎｅｅｒｉｎｇ，Ｖｏｌ．１３４（１９９６），ｐｐ．２２３−２４０．Ｓｕｎ，Ｗ．，Ｃｈａｉｋｏｆ，Ｅ．Ｌ．ａｎｄＬｅｖｅｎｓｔｏｎ，Ｍ．Ｅ．，Ｎｕｍｅｒｉｃａｌａｐｐｒｏｘｉｍａｔｉｏｎｏｆｔａｎｇｅｎｔｍｏｄｕｌｕｓｆｏｒｆｉｎｉｔｅｅｌｅｍｅｎｔｉｍｐｌｅｍｅｎｔａｔｉｏｎｓｏｆｎｏｎｌｉｎｅａｒｈｙｐｅｒｅｌａｓｔｉｃｍａｔｅｒｉａｌｍｏｄｅｌｓ，ＪｏｕｒｎａｌｏｆＢｉｏｍｅｃｈａｎｉｃａｌＥｎｇｉｎｅｅｒｉｎｇ，Ｖｏｌ．１３０，Ｎｏ．６（２００８），ｐｐ．０６１００６３．Ｌａｉ，Ｋ．Ｌ．ａｎｄＣｒａｓｓｉｄｉｓ，J．Ｌ．，Ｅｘｔｅｎｓｉｏｎｓｏｆｔｈｅｆｉｒｓｔａｎｄｓｅｃｏｎｄｃｏｍｐｌｅｘ−ｓｔｅｐｄｅｒｉｖａｔｉｖｅａｐｐｒｏｘｉｍａｔｉｏｎｓ，ＪｏｕｒｎａｌｏｆＣｏｍｐｕｔａｔｉｏｎａｌａｎｄＡｐｐｌｉｅｄＭａｔｈｅｍａｔｅｃｓ，Ｖｏｌ．２１９（２００８），ｐｐ.２７６−２９３．田中真人、藤川正毅、複素数階微分を用いた整合接線剛性の数値近似と大変形問題への応用、日本機械学会論文集Ａ編、Ｖｏｌ．７７,Ｎｏ．７３３（２０１１），ｐｐ．２７−３８．

上記の非特許文献１、２に記載の方法により計算された材料ヤコビアンの近似精度は、数値微分で用いる変動量Δｘの大きさに依存し、大き過ぎると打ち切り誤差が、小さ過ぎると丸め誤差が発生してしまう。変動量Δｘの最適な大きさは、打ち切り誤差と丸め誤差のトレードオフの関係を計りながら決定する必要があるが、材料パラメータ、幾何学的データなどの絶対値に依存して最適なΔｘの値は変化し、明確な指針を得ることは難しい。実際のところ、現状では最適なΔｘの値は経験的にしか評価できていない。それ故、Δｘの値はしばしば“マジックナンバー”と呼ばれる。

また、非特許文献４に記載の複素数階微分法を扱うことができるのは原理的に１階微分のみである。それ以上の高階微分では、前方オイラー法と同じ丸め誤差が生ずる。一般ユーザが新しい構成則を実装する場合を考えると、エネルギー関数式から応力と材料ヤコビアンの両方を同時に求められることが望ましい。すなわち、丸め誤差のない１階および２階の微分近似法と、それを応力および材料ヤコビアンの導出に高効率に適用できる技術が望まれるが、従来技術ではそのような方法は皆無であった。

本発明は、上記の事情を鑑みてなされたもので、誤差の発生を抑制して、関数の１階微分値、２階微分値を計算することができる情報処理装置及びプログラムを提供することを目的とする。

第１の発明に係る情報処理装置は、虚数単位である２つの数ε₁、ε₂であって、ε₁、ε₂のそれぞれを２乗すると０となり、互いに乗算に関して交換可能な数として定義されたε₁、ε₂を用いて、関数に対するテンソルの方向微分を行う情報処理装置であって、入力された２階のテンソル量Ｆの値（Ｆ＝Ｆ＾）に基づいて、以下の（Ｉ）式に従ってΔＦ₁ ^（ij）を計算する第１増分量計算手段と、前記テンソル量Ｆの値（Ｆ＝Ｆ＾）に基づいて、以下の（ＩＩ）式に従って~ΔＦ₂ ^（kl）を計算する第２増分量計算手段と、前記計算されたΔＦ₁ ^（ij）、及び前記計算された~ΔＦ₂ ^（kl）を用いて、関数Ｗ（Ｆ＾＋ΔＦ₁ ^（ij）＋~ΔＦ₂ ^（kl））を計算する関数計算手段と、前記関数計算手段によって計算された関数Ｗ（Ｆ＾＋ΔＦ₁ ^（ij）＋~ΔＦ₂ ^（kl））におけるε₁の係数を取り出して、関数Ｗ（Ｆ）の前記テンソル量Ｆによる１階微分を用いて求められる２階のテンソル量である第１物理量のｉｊ成分とする第１物理量計算手段と、前記関数計算手段によって計算された関数Ｗ（Ｆ＾＋ΔＦ₁ ^（ij）＋~ΔＦ₂ ^（kl））におけるε₁・ε₂の係数を取り出して、前記関数Ｗ（Ｆ）の前記テンソル量Ｆによる２階微分を用いて求められる４階のテンソル量である第２物理量のｉｊｋｌ成分とする第２物理量計算手段と、を含んで構成されている。
ただし、ｅ _ｉ、ｅ _j 、ｅ _k 、ｅ _l はデカルト座標系における単位基底ベクトルであり、
はテンソル積である。

第２の発明に係るプログラムは、虚数単位である２つの数ε₁、ε₂であって、ε₁、ε₂のそれぞれを２乗すると０となり、互いに乗算に関して交換可能な数として定義されたε₁、ε₂を用いて、関数に対するテンソルの方向微分を行うためのプログラムであって、コンピュータを、入力された２階のテンソル量Ｆの値（Ｆ＝Ｆ＾）に基づいて、以下の（ＩＩＩ）式に従ってΔＦ₁ ^（ij）を計算する第１増分量計算手段、前記テンソル量Ｆの値（Ｆ＝Ｆ＾）に基づいて、以下の（ＩＶ）式に従って~ΔＦ₂ ^（kl）を計算する第２増分量計算手段、前記計算されたΔＦ₁ ^（ij）、及び前記計算された~ΔＦ₂ ^（kl）を用いて、関数Ｗ（Ｆ＾＋ΔＦ₁ ^（ij）＋~ΔＦ₂ ^（kl））を計算する関数計算手段、前記関数計算手段によって計算された関数Ｗ（Ｆ＾＋ΔＦ₁ ^（ij）＋~ΔＦ₂ ^（kl））におけるε₁の係数を取り出して、関数Ｗ（Ｆ）の前記テンソル量Ｆによる１階微分を用いて求められる２階のテンソル量である第１物理量のｉｊ成分とする第１物理量計算手段、及び前記関数計算手段によって計算された関数Ｗ（Ｆ＾＋ΔＦ₁ ^（ij）＋~ΔＦ₂ ^（kl））におけるε₁・ε₂の係数を取り出して、前記関数Ｗ（Ｆ）の前記テンソル量Ｆによる２階微分を用いて求められる４階のテンソル量である第２物理量のｉｊｋｌ成分とする第２物理量計算手段として機能させるためのプログラムである。
ただし、ｅ _ｉ、ｅ _j 、ｅ _k 、ｅ _l はデカルト座標系における単位基底ベクトルであり、
はテンソル積である。

第１の発明及び第２の発明によれば、第１増分量計算手段によって、入力された２階のテンソル量Ｆの値（Ｆ＝Ｆ＾）に基づいて、（Ｉ）式に従ってΔＦ₁ ^（ij）を計算する。第２増分量計算手段によって、前記テンソル量Ｆの値（Ｆ＝Ｆ＾）に基づいて、（ＩＩ）式に従って~ΔＦ₂ ^（kl）を計算する。

そして、関数計算手段によって、前記計算されたΔＦ₁ ^（ij）、及び前記計算された~ΔＦ₂ ^（kl）を用いて、関数Ｗ（Ｆ＾＋ΔＦ₁ ^（ij）＋~ΔＦ₂ ^（kl））を計算する。第１物理量計算手段によって、前記関数計算手段によって計算された関数Ｗ（Ｆ＾＋ΔＦ₁ ^（ij）＋~ΔＦ₂ ^（kl））におけるε₁の係数を取り出して、関数Ｗ（Ｆ）の前記テンソル量Ｆによる１階微分を用いて求められる２階のテンソル量である第１物理量のｉｊ成分とする。第２物理量計算手段によって、前記関数計算手段によって計算された関数Ｗ（Ｆ＾＋ΔＦ₁ ^（ij）＋~ΔＦ₂ ^（kl））におけるε₁・ε₂の係数を取り出して、前記関数Ｗ（Ｆ）の前記テンソル量Ｆによる２階微分を用いて求められる４階のテンソル量である第２物理量のｉｊｋｌ成分とする。

このように、関数Ｗ（Ｆ＾＋ΔＦ₁ ^（ij）＋~ΔＦ₂ ^（kl））におけるε₁の係数を取り出して、関数Ｗ（Ｆ）のテンソル量Ｆによる１階微分を用いて求められる第１物理量のｉｊ成分とし、関数Ｗ（Ｆ＾＋ΔＦ₁ ^（ij）＋~ΔＦ₂ ^（kl））におけるε₁・ε₂の係数を取り出して、関数Ｗ（Ｆ）のテンソル量Ｆによる２階微分を用いて求められる第２物理量のｉｊｋｌ成分とすることにより、誤差の発生を抑制して、関数の１階微分を用いて求められる物理量、２階微分を用いて求められる物理量を計算することができる。

上記の第１の発明において、前記関数は、シミュレーション対象に関する関数であって、前記第１物理量計算手段は、シミュレーションに用いられる前記第１物理量を計算し、前記第２物理量計算手段は、シミュレーションに用いられる前記第２物理量を計算するようにすることができる。

また、上記の第１の発明に係る情報処理装置は有限要素法（ＦＥＭ）を用いたシミュレーションを行うシミュレーション手段を更に含み、前記テンソル量Ｆは、ひずみを示す変形勾配テンソルであって、前記シミュレーションは、材料の挙動に関するシミュレーションであって、前記第１物理量は、応力テンソルであり、前記第２物理量は、材料ヤコビアンであり、前記シミュレーション手段は、前記第１物理量計算手段によって計算された応力テンソル、及び前記第２物理量計算手段によって計算された材料ヤコビアンを用いて、シミュレーションを行うようにすることができる。

以上説明したように、本発明の情報処理装置及びプログラムによれば、誤差の発生を抑制して、関数の１階微分値、２階微分値を計算することができる、という効果が得られる。

本発明の第１の参考例に係る情報処理装置を示すブロック図である。本発明の第１の参考例に係る情報処理装置の微分計算処理ルーチンの内容を示すフローチャートである。本発明の第１の実施の形態に係る情報処理装置のシミュレーション処理ルーチンの内容を示すフローチャートである。本発明の第１の実施の形態に係る情報処理装置によるエネルギー関数をＦＥＭ計算に組み込む処理の処理ルーチンの内容を示すフローチャートである。

以下、図面を参照して本発明の実施の形態を詳細に説明する。

図１に示すように、第１の参考例に係る情報処理装置１０は、ＣＰＵ１２、ＲＯＭ１４、ＲＡＭ１６、ＨＤＤ１８、通信インタフェース２０、及びこれらを相互に接続するためのバス２２を備えている。

ＣＰＵ１２は、各種プログラムを実行する。ＲＯＭ１４には、各種プログラムやパラメータ等が記憶されている。ＲＡＭ１６は、ＣＰＵ１２による各種プログラムの実行時におけるワークエリア等として用いられる。記録媒体としてのＨＤＤ１８には、後述する微分計算処理ルーチンを実行するためのプログラムを含む各種プログラムや各種データが記憶されている。

第１の参考例における情報処理装置１０による微分計算方法では、後述するＭｕｌｔｉｄｕａｌｎｕｍｂｅｒを用いて、１変数の関数の１階微分値および２階微分値を計算する。

ここで、Ｍｕｌｔｉｄｕａｌｎｕｍｂｅｒを用いた微分計算の原理について説明する。

まず、Ｍｕｌｔｉｄｕａｌｎｕｍｂｅｒを以下のように定義する。

Ｍｕｌｔｉｄｕａｌｎｕｍｂｅｒとは複素数の変種であり、ε_１、ε_２の２種類の虚数単位を持ち、以下の性質を有するとする。

すなわち、それぞれの虚数単位を二乗すると０になり、２種類の虚数単位は、互いに乗算に関して交換可能である。以上のＭｕｌｔｉｄｕａｌｎｕｍｂｅｒに対して、通常の複素数と同様に、四則演算や初等関数の定義を自然に拡張することができる。以下に、Ｍｕｌｔｉｄｕａｌｎｕｍｂｅｒの代表的な演算例を示す。（ａ_ｉ,ｂ_ｉ（ｉ＝１，２，３，４）は全て実数とする。）

以上の演算規則を持つＭｕｌｔｉｄｕａｌｎｕｍｂｅｒを、以下の（３）式に示す関数ｆ（ｘ）に対するＴａｙｌｏｒ展開式の微小変位量Δｘに代入すると、以下の（４）式が得られる。

すなわち、関数ｆ（ｘ）のｘ＝ａにおける微分値を計算したい場合は、まず、関数ｆ（ｘ）にｘ＝ａ＋ε_１＋ε_２を代入し、関数計算をＭｕｌｔｉｄｕａｌｎｕｍｂｅｒへ機械的に置き換える。その計算結果に対して、ε_１もしくはε_２の係数を取り出すと１階の微分値ｆ’（ａ）が、ε_１ε_２の係数を取り出すと２階の微分値ｆ'’（ａ）が自動的に得られる。

次に、第１の参考例に係る情報処理装置１０による微分計算を行う際の動作について説明する。

まず、情報処理装置１０に、関数ｆ（ｘ）及び微分値を算出するときの変数ｘの値ａが入力されると、情報処理装置１０によって、図２に示す微分計算処理ルーチンが実行される。

まず、ステップ１００において、入力された関数ｆ（ｘ）に、ｘ＝ａ＋ε₁＋ε₂を代入して、関数ｆ（ａ＋ε₁＋ε₂）を計算する。

そして、ステップ１０２で、上記ステップ１００の計算結果から、ε₁又はε₂の係数を取り出して、１階の微分値ｆ’（ａ）を出力する。また、ステップ１０４で、上記ステップ１００の計算結果から、ε₁・ε₂の係数を取り出して、２階の微分値ｆ’’（ａ）を出力し、微分計算処理ルーチンを終了する。

以上説明したように、第１の参考例に係る情報処理装置によれば、Ｍｕｌｔｉｄｕａｌｎｕｍｂｅｒを用いて、関数ｆ（ａ＋ε₁＋ε₂）を計算し、関数ｆ（ａ＋ε₁＋ε₂）におけるε₁またはε₂の係数を、関数をスカラー量で微分したときの１階微分値ｆ’（ａ）として取り出すと共に、関数ｆ（ａ＋ε₁＋ε₂）におけるε₁・ε₂の係数を、２階微分値ｆ’’（ａ）として取り出すことにより、誤差の発生を抑制して、関数の１階微分値、２階微分値を計算することができる。

なお、上記の参考例では、四則演算の定義を説明したが、Ｍｕｌｔｉｄｕａｌｎｕｍｂｅｒを引数にとった初等関数の計算も可能となる。以下にいくつかのＭｕｌｔｉｄｕａｌｎｕｍｂｅｒによる初等関数の計算例を示す。

次に、第２の参考例について説明する。なお、第２の参考例に係る情報処理装置は、第１の参考例と同様の構成となるため、同一符号を付して説明を省略する。

第２の参考例では、２変数の関数の偏微分値を計算している点が、第１の参考例と異なっている。

２変数の関数の偏微分値は、１変数の関数と同様に、以下のように自然に拡張できる。すなわち、以下の（５）式に示す２変数の関数に対して、微小変位量Δｘにε_１を代入し、微小変位量Δｙにε_２を代入すると、以下の（６）式が得られる。

ここで、関数ｇ（ｘ、ｙ）のｘ＝ａ、ｙ＝ｂにおける微分値を計算したい場合は、関数ｇ（ｘ、ｙ）もｘ＝ａ＋ε_１、ｙ＝ｂ＋ε_２を代入し、計算された結果に対して、ε_１の係数を取り出すと１階の偏微分値∂ｇ（ａ、ｂ）／∂ｘが、ε_２の係数を取り出すと１階の偏微分値∂ｇ（ａ、ｂ）／∂ｙが、ε_１ε_２の係数を取り出すと２階の偏微分値∂^２ｇ（ａ、ｂ）／∂ｘ∂ｙが自動的に得られる。

次に、第２の参考例に係る情報処理装置１０による微分計算を行う際の動作について説明する。

まず、情報処理装置１０に、関数ｇ（ｘ、ｙ）及び偏微分値を算出するときの変数ｘ，ｙの値ａ、ｂが入力されると、情報処理装置１０によって、上記図２と同様の微分計算処理ルーチンが実行される。

まず、入力された関数ｇ（ｘ、ｙ）に、ｘ＝ａ＋ε₁、ｙ＝ｂ＋ε₂を代入して、関数ｇ（ａ＋ε₁、ｂ＋ε₂）を計算する。

そして、上記の計算結果から、ε₁の係数を取り出して、１階の偏微分値∂ｇ（ａ、ｂ）／∂ｘを出力する。また、上記の計算結果から、ε₂の係数を取り出して、１階の偏微分値∂ｇ（ａ、ｂ）／∂ｙを出力する。また、上記の計算結果から、ε₁・ε₂の係数を取り出して、２階の偏微分値∂²ｇ（ａ、ｂ）／∂ｘ∂ｙを出力し、微分計算処理ルーチンを終了する。

以上説明したように、第２の参考例に係る情報処理装置によれば、Ｍｕｌｔｉｄｕａｌｎｕｍｂｅｒを用いて、ｇ（ａ＋ε₁、ｂ＋ε₂）を計算し、関数ｇ（ａ＋ε₁、ｂ＋ε₂）におけるε₁の係数を、１階偏微分値∂ｇ（ａ、ｂ）／∂ｘとして取り出し、関数ｇ（ａ＋ε₁、ｂ＋ε₂）におけるε₂の係数を、１階偏微分値∂ｇ（ａ、ｂ）／∂ｙとして取り出し、関数ｇ（ａ＋ε₁、ｂ＋ε₂）におけるε₁・ε₂の係数を、２階偏微分値∂²ｇ（ａ、ｂ）／∂ｘ∂ｙとして取り出すことにより、誤差の発生を抑制して、関数の１階偏微分値、２階偏微分値を計算することができる。

次に、第１の実施の形態について説明する。なお、第１の実施の形態に係る情報処理装置は、第１の参考例と同様の構成となるため、同一符号を付して説明を省略する。

第１の実施の形態では、ＦＥＭを用いたシミュレーションを行っている点と、テンソル方向微分による微分値を計算している点が、第１の参考例と異なっている。

第１の実施の形態における情報処理装置１０による微分計算方法では、Ｍｕｌｔｉｄｕａｌｎｕｍｂｅｒを用いて、エネルギー関数のテンソル方向微分による１階微分値および２階微分値を計算する。また、情報処理装置１０による材料のシミュレーション方法では、上記の微分計算方法により計算される１階微分値および２階微分値を用いて、ＦＥＭ計算を行い、入力されるひずみ（テンソル量）に対する応力を、シミュレーション結果として計算する。

次に、Ｍｕｌｔｉｄｕａｌｎｕｍｂｅｒを用いた応力および材料ヤコビアンの自動計算の原理について説明する。

以下にエネルギー関数から応力と材料ヤコビアンを計算する方法を示す。通常、汎用ソフトウェアをはじめとしたＦＥＭプログラムにおける材料構成則のユーザサブルーチンでは、「参照のため引き渡される変数」として、変形勾配テンソルＦが入力される。ユーザはＦを用いて、エネルギー関数から計算されるＣａｕｃｈｙ応力σと材料ヤコビアンＣ^∇MJ（４階テンソル）の各成分を引き渡すプログラムを実装することになる。また、汎用ＦＥＭソフトウェアの材料構成則ユーザサブルーチンでは、材料ヤコビアンにＫｉｒｃｈｈｏｆｆ応力τのＪａｕｍａｎｎ速度を用いて、ｕｐｄａｔｅｄＬａｇｒａｎｇｅ法による定式化を採用していることが多い。ここでも、その定式化をもとに説明する。Ｋｉｒｃｈｈｏｆｆ応力τ、τのＪａｕｍａｎｎ速度τ^∇Ｊ、対応する材料ヤコビアンＣ^∇ＭＪの定義式はそれぞれ以下の（７）式〜（９）式となる。

ここで、「^・」は物質時間微分、「：」はテンソルの二組の基底ベクトルに対する縮約を表す。また、Ｊは体積変化率であり、変形勾配テンソルＦを用いて、以下の（１０）式のように表される。

また、Ｄ、Ｗは、以下の（１１）式の変形速度勾配テンソルＬの対称成分、反対称成分である。

ここで、Ｔ^−１はテンソルＴの逆行列を表す。

次に、応力の計算法について説明する。

τの対称性を考慮すれば、τのｉｊ成分τ_ｉｊは、以下の（１２）式のように表すことができる。

ここで、ｅ_ｉはデカルト座標系における単位基底ベクトルであり、
はテンソル積である。まず、Ｗ（Ｆ）のＦによる微分（テンソルの方向微分）を考える。のちのτの導出過程を単純化するために、変形勾配テンソルＦの微小増分量をΔＦ_１ ^（ｉｊ）とすれば、以下の（１３）式に示す近似式が得られる。

いま微小増分量ΔＦ_１ ^（ｉｊ）を、Ｍｕｌｔｉｄｕａｌｎｕｍｂｅｒの虚数単位ε_１を用いて以下の（１４）式のように定義する。

次に、上記（１４）式を上記（１３）式に代入して右辺を整理すれば、以下の（１５）式を得る。

ただし、Ｔ^ＴはテンソルＴの転置を表す。上記（１５）式の右辺には、以下の（１６）式に示す第１Ｐｉｏｌａ−Ｋｉｒｃｈｈｏｆｆ応力Ｐが含まれており、τとの関係は以下の（１６）式のようである。

なお、応力τが、関数Ｗ（Ｆ）のテンソル量Ｆによる１階微分に基づく第１物理量の一例である。また、材料ヤコビアンが、関数Ｗ（Ｆ）のテンソル量Ｆによる２階微分に基づく第２物理量の一例である。

上記（１２）式、（１６）式を用いて、上記（１５）式を整理すれば、Ｗ（Ｆ）からτを算出する、以下の（１７）式が得られる。

ここで、
はε_１の係数を取り出す演算子とする。

次に、材料ヤコビアンの計算法について説明する。

エネルギー関数Ｗ（Ｆ）から材料ヤコビアンＣ^∇ＭＪを計算する方法を示す。Ｋｉｒｃｈｈｏｆｆ応力τから、材料ヤコビアンＣ^∇ＭＪを計算する方法は、以下の（１８）式の通りである。

ここで、
はε_２の係数を取り出す演算子であり、ΔＦ_２ ^（ｋｌ）は以下の（１９）式のように定義される。

上記（１７）式と（１８）式とを組み合わせると、以下の（２０）式が得られる。

ここで、

はε_１ε_２の係数を取り出す演算子であり、~ΔＦ_２ ^（ｋｌ）は以下の（２１）式のように定義される。

このとき、応力は以下の（２２）式で求まる。

また、増分量ΔＦ_１ ^（ｉｊ）、~ΔＦ_２ ^（ｋｌ）の求め方について説明する。

以下に示すように、ΔＦ_１ ^（ｉｊ）を、Cauchy応力テンソルσを導くように設定し、~ΔＦ_２ ^（ｋｌ）を、材料ヤコビアンＣ^∇MJを導くように設定している．

Cauchy応力テンソルσとKirchhoff応力テンソルτは、上記（７）式で関係付けられている。すなわち、Kirchhoff応力テンソルτが求められれば、Ｊで割ることでただちにCauchy応力テンソルσを求めることができる。したがって、エネルギー関数Ｗ（Ｆ）から Kirchhoff応力テンソルτを導くΔＦ_１ ^（ｉｊ）の設定方法を以下に示す．

上記（１６）式で示しているように，エネルギー関数Ｗ（Ｆ）とKirchhoff応力テンソルτには、以下の（２３）式に示す関係がある。

なお、上記（２３）式が、第１物理量と関数Ｗ（Ｘ）との関係式に対応している。

τの対称性を考えて、上記（２３）式の両辺の転置をとると、以下の（２４）式が得られる。

上記（１２）式の方法で、上記（２４）式のτのｉｊ成分τ_ijを求めると、以下の（２５）式となる。

さらに、上記（２５）式を、以下の（２６）式のように式変形する。

ここで、テンソルの方向微分とＭｕｌｔｉｄｕａｌｎｕｍｂｅｒによる微分法の関係を以下に示す。まず、スカラー値テンソル関数に対するテンソルの方向微分の定義を示す。テンソルを独立変数とするスカラー関数を「スカラー値テンソル関数」と呼ぶ。いま、スカラー値テンソル関数Ｇ（Ａ）を考え（ここではＧはスカラー、Ａは2階のテンソル）、このＧをΔＡの方向にＡで微分することを以下の（２７）式で表す。

ここで、記号ＤＧ（Ａ）［ΔＡ］で方向微分を表し、ΔＡは方向テンソルと呼ばれる。上記（２７）式の第2辺の微分の式を、Ｍｕｌｔｉｄｕａｌｎｕｍｂｅｒを用いて書き直すと、以下の（２８）式のようになる。

上記（２７）式、（２８）式より、以下の（２９）式となる。

上記（２７）式〜（２９）式が、テンソル方向微分とε₁による微分との関係に対応している。

上記（２９）式と上記（２６）式とを比べると、ＧにＷを、ＡにＦを、ΔＡに

をそれぞれ代入すると、τ_ijは以下の（３０）式で表わされることがわかる。

すなわち、

とおくと、上記（１７）式によってτ_ijが求まることがわかる。

次に、材料ヤコビアンＣ^▽ＭＪを導くための~ΔＦ_２ ^（ｋｌ）の設定方法を以下に示す。上記（９）式のようにＣ^▽ＭＪはτ^▽JとＤの間の係数として定義される。まず、上記（８）式、（９）式を、以下の（３１）式、（３２）式のように増分形にする。

ここで、ΔＤ、ΔＷは変形勾配テンソルの増分形ΔＦを用いて、以下の（３３）式、（３４）式のように表される。

いま、ΔＦに

を代入すると、上記（３１）式、（３２）式より、以下の（３５）式〜（３８）式が導かれる。

なお、上記（３７）式が、第１物理量の増分と第２物理量との関係式に対応している。

なお、上記（３９）式が、テンソル方向微分とε₂による微分との関係に対応している。

Ｃ^▽ＭＪの対称性より、上記（３９）式を成分表示すると、以下の（４０）式となる。

上記（４０）式の左辺の微分の式を、Ｍｕｌｔｉｄｕａｌｎｕｍｂｅｒを用いて書き直すと以下の（４１）式のようになる。

なお、上記（４１）式が、第２物理量と前記第１物理量との関係式に対応している。

したがって、

とおくと、上記（１８）式で、Ｃ^▽ＭＪのｉｊｋｌ成分（Ｃ^▽ＭＪ）_ijklが得られることがわかる。

最後に、上記（３０）式と（４１）式とを同時に評価できるように式をさらに工夫する。なお、上記（３０）式が、第1物理量と関数Ｗ（Ｘ）との関係式に対応している。

ΔＦ_１ ^（ｉｊ）は応力の計算(エネルギー関数Ｗを1階微分)のために使われ、ΔＦ_２ ^（ｋｌ）は材料ヤコビアンの計算(エネルギー関数Ｗを2階微分)のために使われていることを考えると、ε₁の係数に応力の成分を、ε₁ε₂の係数に材料ヤコビアンの成分が現れるように設定するのが自然である。そのため、ΔＦ_２ ^（ｋｌ）の代わりに、上記（２１）式の~ΔＦ_２ ^（ｋｌ）を使用すると、この目的がうまく果たされる。実際に上記（２０）式、（２２）式を計算して確かめてみる。

上記（４２）式のように、上記（２０）式が導出され、上記（４３）式のように、上記（２２）式が導出された。なお、上記（４２）式が、第２物理量と関数Ｗ（Ｘ）との関係式に対応している。

ここで、ＦＥＭ計算の概要について説明する。有限要素法は、構造解析などにおいて、変形に対して無限の自由度を持つ物体を有限の自由度を持つ有限要素の集合体、すなわち、微小な部分の集合体として近似し、集合体に対して成立する連立一次方程式を解く方法である。この微小な部分を有限要素と呼ぶ。有限要素は節点と呼ばれる点を結合したものとして規定される。ある要素は他の要素と節点により結合される。力も節点を通じて結合されている。いかに複雑な形状の構造物であっても、これらを有限要素に区分けし、有限要素の集合体として表すことができる。各有限要素は、材料の挙動を表す剛性マトリクス（非線形解析であれば接線剛性マトリクス）を有する。

次に、第１の実施の形態に係る情報処理装置１０によるＦＥＭを用いたシミュレーションを行う際の動作について説明する。

まず、情報処理装置１０に、材料のシミュレーションで用いる、材料科学の分野で提案された複数種類のエネルギー関数Ｗ（Ｆ）が入力される。また、エネルギー関数Ｗ（Ｆ）から、Cauchy応力テンソルσを導くためのΔＦ_１ ^（ｉｊ）の式と、材料ヤコビアンＣ^∇MJを導くための~ΔＦ_２ ^（ｋｌ）の式を、予め計算しておき、情報処理装置１０に入力する。そして、情報処理装置１０によって、図３に示すシミュレーション処理ルーチンが実行される。

まず、ステップ３００において、情報処理装置１０に、材料の応力ひずみ曲線の実験データが入力されたか否かを判定する。当該実験データが入力されると、ステップ３０２へ進み、入力されたエネルギー関数のうちの何れか１つを、シミュレーションで用いるエネルギー関数として設定する。また、当該エネルギー関数に含まれる材料パラメータを設定する。たとえば、当該エネルギー関数に含まれる材料パラメータを、上記ステップ３００で入力された応力ひずみ曲線の実験データと合致するように同定する。

そして、ステップ３０４において、上記ステップ３０２で設定されたエネルギー関数をＦＥＭ計算に組み込む処理を行う。

上記ステップ３０４は、図４に示す処理ルーチンによって実現される。

ステップ３１０において、情報処理装置１０に、ひずみのテンソル量（変形勾配テンソル）Ｆ＾が入力されたか否かを判定する。そして、ステップ３１２において、エネルギー関数ＷからCauchy応力を計算するように予め設定され、かつ、入力されたΔＦ_１ ^（ｉｊ）の式を、上記ステップ３１０で入力された変形勾配Ｆ＾に基づいて、（ｉｊ）成分毎に計算する。このとき、ＭＤＮ（ＭｕｌｔｉｄｕａｌＮｕｍｂｅｒｓ）であるε₁を用いて上記（１４）式のようにΔＦ_１ ^（ｉｊ）が計算される。

次のステップ３１４では、エネルギー関数Ｗから材料ヤコビアンＣ^∇MJを計算するように予め設定され、かつ、入力された~ΔＦ_２ ^（ｋｌ）の式を、上記ステップ３１０で入力された変形勾配Ｆ＾に基づいて、（ｋｌ）成分毎に計算する。このとき、ＭＤＮであるε₁、ε₂を用いて上記（２１）式のように~ΔＦ_２ ^（ｋｌ）が計算される。

そして、ステップ３１６において、上記ステップ３１２の計算結果及び上記ステップ３１４の計算結果に基づいて、（ｉｊ）成分，（ｋｌ）成分の組み合わせ毎に、エネルギー関数Ｗ（Ｆ＾＋ΔＦ_１ ^（ｉｊ）＋~ΔＦ_２ ^（ｋｌ））の計算を行う。

次のステップ３１８では、後述するＦＥＭ計算におけるCauchy応力の成分σ_ijの計算を、上記ステップ３１６で計算したエネルギー関数Ｗ（Ｆ＾＋ΔＦ_１ ^（ｉｊ）＋~ΔＦ_２ ^（ｋｌ））のε₁の係数を取り出す処理に置き換える。

また、ステップ３２０において、後述するＦＥＭ計算における材料ヤコビアンの成分（Ｃ^∇MJ）_ijklの計算を、上記ステップ３１６で計算したエネルギー関数Ｗ（Ｆ＾＋ΔＦ_１ ^（ｉｊ）＋~ΔＦ_２ ^（ｋｌ））のε₁ε₂の係数を取り出す処理に置き換えて、当該処理ルーチンを終了する。

これによって、エネルギー関数式から応力と材料ヤコビアンを自動的に算出する処理を、FEMプログラムの中の材料構成則のルーチンの中へ組み込むことができる。

そして、シミュレーション処理ルーチンのステップ３０６において、ＦＥＭ計算により、積分点ごとに、上記ステップ３１０で入力された変形勾配テンソルＦに対する応力を計算する。このとき、応力と材料ヤコビアンの計算では、上記ステップ３１８、３２０で置換された計算方法を用いる。

ここで、ＦＥＭ計算の流れについて説明する。なお、以下では、荷重制御の増分反復（インクリメント−イタレーション）型の非線形ＦＥＭの流れについて説明する。

（１）解析対象を有限要素に区分けする。
（２）境界条件を定める。
（３）各有限要素の要素接線剛性マトリクスを計算する。
（４）要素接線剛性マトリクスを重ねて併せて全体接線剛性マトリクスを計算する。
（５）拘束した変位自由度に対応する全体接線剛性マトリクスの成分を消去し、全体接線剛性マトリクスを縮退する。
（６）荷重増分を与える。
（７）連立一次方程式を解いて、変位増分を計算する。
（８）上記（７）で計算した変位増分を全変位量に加算し、全変位量を更新する。
（９）全変位量から各有限要素のひずみ、応力を計算する。
（１０）各有限要素の応力から、各有限要素の等価節点力を計算する。
（１１）各有限要素の等価節点力を重ね合わせて、構造全体の等価節点力を計算する。
（１２）上記（１１）で計算された構造全体の等価節点力と上記（６）で与えた荷重増分を比べて力が釣り合っているか確認する。
（１３）力が釣り合っていなければ、上記（３）に戻り、上記（７）で計算した変位増分を加味した要素接線剛性マトリクスを計算する。
（１４）力が釣り合うまで、上記（３）〜（１３）の計算ステップを繰り返す（この反復計算をＮｅｗｔｏｎ−Ｒａｐｈｓｏｎ反復（イタレーション）法と呼ぶ）。
（１５）力が釣り合ったら、次の荷重増分を加え、上記（３）〜（１５）の計算を繰り返す（この増分計算をインクリメントと呼ぶ）。
（１６）荷重を増分し続け、所望の荷重値になったら計算を終了する。
（１７）ポスト処理により、全体の変位分布、荷重分布、ひずみ分布、応力分布を表示する。

また、以下では、上記（３）における要素接線剛性マトリクスの作成方法について示す。

要素接線剛性マトリクスの算出に際しては、有限要素の体積内での積分が必要となり、通常これには数値積分（Ｇａｕｓｓ積分、Ｎｅｗｔｏｎ−Ｃｏｔｅｓ積分など）を用いる。すなわち、要素内の複数の積分点において剛性マトリクスを算出し、それに重みをかけて総和を取る。さらにアイソパラメトリック四辺形要素など、一般座標から自然座標への写像を必要とする有限要素の場合、各積分点におけるヤコビ行列（ヤコビアン）の算出を行う。

（３−１）積分点座標と重みを求める。
（３−２）積分点ごとに、ヤコビアンとその逆行列を求める。
（３−３）積分点ごとに、変位ひずみマトリクスを求める。
（３−４）積分点ごとに、応力、ひずみを求める。
（３−５）積分点ごとに、材料ヤコビアンを求める。
（３−６）積分点ごとに、接線剛性マトリクスを求める。
（３−７）各積分点の接線剛性マトリクスに重みを掛けて総和をとり、要素接線剛性マトリクスを計算する。

上記ステップ３２０の置換では、上記（３−４）の応力の計算、上記（３−５）の材料ヤコビアンの計算における計算方法が置換される。

そして、ステップ３０８において、上記ステップ３０６で計算された応力と、上記ステップ３００で入力された実験データにおける、上記ステップ３１０で入力されたひずみ量Ｆ＾に対する応力とを比較し、実験値とＦＥＭの計算値とが一致しているか否かを判定する。実験値とＦＥＭの計算値とが一致している場合には、このときのエネルギー関数を出力して、シミュレーション処理ルーチンを終了する。一方、実験値とＦＥＭの計算値とが一致していない場合には、上記ステップ３０２へ戻り、別のエネルギー関数を、シミュレーションで用いるエネルギー関数として設定する。

以上説明したように、第１の実施の形態に係る情報処理装置は、Ｍｕｌｔｉｄｕａｌｎｕｍｂｅｒを用いて、エネルギー関数Ｗ（Ｆ＾＋ΔＦ₁ ^（ij）＋~ΔＦ₂ ^（kl））を計算し、エネルギー関数Ｗ（Ｆ＾＋ΔＦ₁ ^（ij）＋~ΔＦ₂ ^（kl））におけるε₁の係数を取り出して、エネルギー関数Ｗ（Ｆ）のテンソル量Ｆによる１階微分に基づく応力を計算し、エネルギー関数Ｗ（Ｆ＾＋ΔＦ₁ ^（ij）＋~ΔＦ₂ ^（kl））におけるε₁・ε₂の係数を取り出して、エネルギー関数Ｗ（Ｆ）のテンソル量Ｆによる２階微分に基づく材料ヤコビアンを計算することにより、誤差の発生を抑制して、応力及び材料ヤコビアンを計算することができる。

Ｍｕｌｔｉｄｕａｌｎｕｍｂｅｒという新しい数の体系に、関数の1階微分値と2階微分値を自動的に計算させる性質があり、さらに、Ｍｕｌｔｉｄｕａｌｎｕｍｂｅｒは、テンソルの方向微分と相性が良く、これを組み合わせることで、エネルギー関数式から、応力と材料ヤコビアンの両方を自動的に計算させることが可能となる。従来、これらの計算は全て手計算によって解析的に導出せねばならなかった。そのため、複雑な材料構成則を扱う場合は、その計算を遂行するための深い専門知識と多大な工程が必要とされた。本実施の形態に係る情報処理装置の導入により、誰でも簡単に短時間でその材料構成則を正確に実装することができる。これにより、ユーザサブルーチン等を通して、材料科学の分野で提案されるエネルギー関数をその複雑さに関係なく容易に汎用有限要素法ソフトウェア上に実装することが可能になり、材料開発のスピードが大幅に向上される。

また、複雑なエネルギー関数式をFEM計算へ容易に組込むことが可能となり、材料試験結果との比較や、材料のまだ観察されていない任意の変形挙動を予測できるようになる。

ゴムなどの高分子エラストマーに代表される超弾性材料は、ひずみエネルギー密度関数Ｗを有する。ひずみエネルギー密度関数Ｗは、物体の変形により蓄積される単位体積当たりの弾性エネルギーとして定義され、等温状態におけるＷの変化は、系の自由エネルギー変化と等価である。したがって、Ｗの関数形が既知であれば任意の変形状態に対する応力−ひずみ関係を求めることが可能となる。複雑な変形下におけるエラストマーの力学応答を知るためにＦＥＭを使用する場合、解析結果の信頼性はＷの関数形に大きく依存する。高分子物理の分野では、分子論的な考察によりエラストマーのネットワーク構造を反映したひずみエネルギー密度関数Ｗが多数提案されている。このようなひずみエネルギー密度関数Ｗを用いると、実験が簡便な一軸引張試験のみの実験データを用いて、多軸変形場の力学応答を高精度に予測することが可能になる。ただし、高精度なひずみエネルギー密度関数Ｗの関数形は複雑化する傾向にあり、これらのＷから応力テンソル、材料ヤコビアンを計算することがＦＥＭへの実装上の障壁となっている。本実施の形態に係る情報処理装置を用いれば、Ｗから応力テンソル、材料ヤコビアンを自動的に算出させることが可能になる。

なお、上記の第１の実施の形態では、材料シミュレーションで用いるエネルギー関数を入力する場合を例に説明したが、これに限定されるものではなく、他のシミュレーションで用いる他の関数を入力するようにしてもよい。たとえば、高分子、金属、非鉄金属、半導体、セラミックス、土質、レオロジー物質、圧電材料、磁性材料、超電導物質、またはこれらを組み合わせた複合材料に対して、様々な応力場における変形量のシミュレーションを行うようにし、当該シミュレーションで用いる関数を入力するようにしてもよい。

本発明のプログラムは、記憶媒体に格納して提供するようにしてもよい。

１０情報処理装置
１２ＣＰＵ
１４ＲＯＭ
１６ＲＡＭ
１８ＨＤＤ

Claims

虚数単位である２つの数ε₁、ε₂であって、ε₁、ε₂のそれぞれを２乗すると０となり、互いに乗算に関して交換可能な数として定義されたε₁、ε₂を用いて、関数に対するテンソルの方向微分を行う情報処理装置であって、
入力された２階のテンソル量Ｆの値（Ｆ＝Ｆ＾）に基づいて、以下の（Ｉ）式に従ってΔＦ₁ ^（ij）を計算する第１増分量計算手段と、
前記テンソル量Ｆの値（Ｆ＝Ｆ＾）に基づいて、以下の（ＩＩ）式に従って~ΔＦ₂ ^（kl）を計算する第２増分量計算手段と、
前記計算されたΔＦ₁ ^（ij）、及び前記計算された~ΔＦ₂ ^（kl）を用いて、関数Ｗ（Ｆ＾＋ΔＦ₁ ^（ij）＋~ΔＦ₂ ^（kl））を計算する関数計算手段と、
前記関数計算手段によって計算された関数Ｗ（Ｆ＾＋ΔＦ₁ ^（ij）＋~ΔＦ₂ ^（kl））におけるε₁の係数を取り出して、関数Ｗ（Ｆ）の前記テンソル量Ｆによる１階微分を用いて求められる２階のテンソル量である第１物理量のｉｊ成分とする第１物理量計算手段と、
前記関数計算手段によって計算された関数Ｗ（Ｆ＾＋ΔＦ₁ ^（ij）＋~ΔＦ₂ ^（kl））におけるε₁・ε₂の係数を取り出して、前記関数Ｗ（Ｆ）の前記テンソル量Ｆによる２階微分を用いて求められる４階のテンソル量である第２物理量のｉｊｋｌ成分とする第２物理量計算手段と、を含む
情報処理装置。
ただし、ｅ _ｉ、ｅ _j 、ｅ _k 、ｅ _l はデカルト座標系における単位基底ベクトルであり、
はテンソル積である。
前記関数は、シミュレーション対象に関する関数であって、
前記第１物理量計算手段は、シミュレーションに用いられる前記第１物理量を計算し、
前記第２物理量計算手段は、シミュレーションに用いられる前記第２物理量を計算する請求項１記載の情報処理装置。
有限要素法（ＦＥＭ）を用いたシミュレーションを行うシミュレーション手段を更に含み、
前記テンソル量Ｆは、ひずみを示す変形勾配テンソルであって、
前記シミュレーションは、材料の挙動に関するシミュレーションであって、
前記第１物理量は、応力テンソルであり、
前記第２物理量は、材料ヤコビアンであり、
前記シミュレーション手段は、前記第１物理量計算手段によって計算された応力テンソル、及び前記第２物理量計算手段によって計算された材料ヤコビアンを用いて、シミュレーションを行う請求項２記載の情報処理装置。
虚数単位である２つの数ε₁、ε₂であって、ε₁、ε₂のそれぞれを２乗すると０となり、互いに乗算に関して交換可能な数として定義されたε₁、ε₂を用いて、関数に対するテンソルの方向微分を行うためのプログラムであって、
コンピュータを、
入力された２階のテンソル量Ｆの値（Ｆ＝Ｆ＾）に基づいて、以下の（ＩＩＩ）式に従ってΔＦ₁ ^（ij）を計算する第１増分量計算手段、
前記テンソル量Ｆの値（Ｆ＝Ｆ＾）に基づいて、以下の（ＩＶ）式に従って~ΔＦ₂ ^（kl）を計算する第２増分量計算手段、
前記計算されたΔＦ₁ ^（ij）、及び前記計算された~ΔＦ₂ ^（kl）を用いて、関数Ｗ（Ｆ＾＋ΔＦ₁ ^（ij）＋~ΔＦ₂ ^（kl））を計算する関数計算手段、
前記関数計算手段によって計算された関数Ｗ（Ｆ＾＋ΔＦ₁ ^（ij）＋~ΔＦ₂ ^（kl））におけるε₁の係数を取り出して、関数Ｗ（Ｆ）の前記テンソル量Ｆによる１階微分を用いて求められる２階のテンソル量である第１物理量のｉｊ成分とする第１物理量計算手段、及び
前記関数計算手段によって計算された関数Ｗ（Ｆ＾＋ΔＦ₁ ^（ij）＋~ΔＦ₂ ^（kl））におけるε₁・ε₂の係数を取り出して、前記関数Ｗ（Ｆ）の前記テンソル量Ｆによる２階微分を用いて求められる４階のテンソル量である第２物理量のｉｊｋｌ成分とする第２物理量計算手段
として機能させるためのプログラム。
ただし、ｅ _ｉ、ｅ _j 、ｅ _k 、ｅ _l はデカルト座標系における単位基底ベクトルであり、
はテンソル積である。
虚数単位である２つの数ε₁、ε₂であって、ε₁、ε₂のそれぞれを２乗すると０となり、互いに乗算に関して交換可能な数として定義されたε₁、ε₂を用いて、シミュレーション対象である材料に関する関数に対するテンソルの方向微分を行う情報処理装置であって、
ひずみを示す変形勾配テンソルとして入力された２階のテンソル量Ｆの値（Ｆ＝Ｆ＾）に基づいて、以下の（Ｖ）式に従ってΔＦ₁ ^（ij）を計算する第１増分量計算手段と、
前記テンソル量Ｆの値（Ｆ＝Ｆ＾）に基づいて、以下の（ＶＩ）式に従って~ΔＦ₂ ^（kl）を計算する第２増分量計算手段と、
前記計算されたΔＦ₁ ^（ij）、及び前記計算された~ΔＦ₂ ^（kl）を用いて、関数Ｗ（Ｆ＾＋ΔＦ₁ ^（ij）＋~ΔＦ₂ ^（kl））を計算する関数計算手段と、
前記関数計算手段によって計算された関数Ｗ（Ｆ＾＋ΔＦ₁ ^（ij）＋~ΔＦ₂ ^（kl））におけるε₁の係数を取り出して、関数Ｗ（Ｆ）の前記テンソル量Ｆによる１階微分を用いて求められる２階のテンソル量である応力テンソルのｉｊ成分とする第１物理量計算手段と、
前記関数計算手段によって計算された関数Ｗ（Ｆ＾＋ΔＦ₁ ^（ij）＋~ΔＦ₂ ^（kl））におけるε₁・ε₂の係数を取り出して、前記関数Ｗ（Ｆ）の前記テンソル量Ｆによる２階微分を用いて求められる４階のテンソル量である材料ヤコビアンのｉｊｋｌ成分とする第２物理量計算手段と、
前記第１物理量計算手段によって計算された応力テンソル、及び前記第２物理量計算手段によって計算された材料ヤコビアンを用いて、有限要素法（ＦＥＭ）を用いた、前記材料の挙動に関するシミュレーションを行うシミュレーション手段と、
を含む情報処理装置。
ただし、ｅ _ｉ、ｅ _j 、ｅ _k 、ｅ _l はデカルト座標系における単位基底ベクトルであり、
はテンソル積である。
虚数単位である２つの数ε₁、ε₂であって、ε₁、ε₂のそれぞれを２乗すると０となり、互いに乗算に関して交換可能な数として定義されたε₁、ε₂を用いて、シミュレーション対象である材料に関する関数に対するテンソルの方向微分を行うためのプログラムであって、
コンピュータを、
ひずみを示す変形勾配テンソルとして入力された２階のテンソル量Ｆの値（Ｆ＝Ｆ＾）に基づいて、以下の（ＶＩＩ）式に従ってΔＦ₁ ^（ij）を計算する第１増分量計算手段、
前記テンソル量Ｆの値（Ｆ＝Ｆ＾）に基づいて、以下の（ＶＩＩＩ）式に従って~ΔＦ₂ ^（kl）を計算する第２増分量計算手段、
前記計算されたΔＦ₁ ^（ij）、及び前記計算された~ΔＦ₂ ^（kl）を用いて、関数Ｗ（Ｆ＾＋ΔＦ₁ ^（ij）＋~ΔＦ₂ ^（kl））を計算する関数計算手段、
前記関数計算手段によって計算された関数Ｗ（Ｆ＾＋ΔＦ₁ ^（ij）＋~ΔＦ₂ ^（kl））におけるε₁の係数を取り出して、関数Ｗ（Ｆ）の前記テンソル量Ｆによる１階微分を用いて求められる２階のテンソル量である応力テンソルのｉｊ成分とする第１物理量計算手段、
前記関数計算手段によって計算された関数Ｗ（Ｆ＾＋ΔＦ₁ ^（ij）＋~ΔＦ₂ ^（kl））におけるε₁・ε₂の係数を取り出して、前記関数Ｗ（Ｆ）の前記テンソル量Ｆによる２階微分を用いて求められる４階のテンソル量である材料ヤコビアンのｉｊｋｌ成分とする第２物理量計算手段、及び
前記第１物理量計算手段によって計算された応力テンソル、及び前記第２物理量計算手段によって計算された材料ヤコビアンを用いて、有限要素法（ＦＥＭ）を用いた、前記材料の挙動に関するシミュレーションを行うシミュレーション手段
として機能させるためのプログラム。
ただし、ｅ _ｉ、ｅ _j 、ｅ _k 、ｅ _l はデカルト座標系における単位基底ベクトルであり、
はテンソル積である。