(57)【要約】【課題】複合材料の材料定数の算出を、効率よく短時間に算出する。また、複合材料中の材料の体積分率の算出を効率よく短時間に行う。【解決手段】複合材料について、母相中に、複数の種類の材料が既知の体積分率で球状粒子として分散した仮想複合材料を定めることにより、この仮想複合材料の有効材料定数を未知数とする非線形方程式を用意する。次に、この方程式を解いて、仮想複合材料の有効材料定数を、複合材料の有効材料定数として求める。このとき、方程式は、セルフコンシステント法を用いて得られる再帰型の式である。この方程式を用いて複合材料中の材料の体積分率を算出する。【選択図】  図１
2010205254.tif 000002
【特許請求の範囲】
【請求項１】
力学的材料定数が既知の複数の種類の材料を、力学的材料定数が既知の母相中に含む複合材料に関して、複合材料の力学的有効材料定数をコンピュータが算出する方法であって、
前記母相中に、前記複数の種類の材料のそれぞれが既知の体積分率で球状粒子として分散した仮想複合材料を、前記複合材料として定めることにより、この仮想複合材料の有効材料定数を未知数とする非線形方程式を用意するステップと、
前記用意した非線形方程式を解いて、前記仮想複合材料の有効材料定数を、前記複合材料の有効材料定数として求めるステップと、を有し、
前記非線形方程式は、前記仮想複合材料中の前記複数の種類の材料のそれぞれに対応する前記球状粒子を囲む周囲の材料定数を、算出すべき前記仮想複合材料の有効材料定数として定めることにより得られる再帰型非線型方程式である複合材料の有効材料定数の算出方法。
【請求項２】
前記母相の材料定数をＣAとし、前記複数の種類の材料の材料定数をＣB1，ＣB2，………ＣBn(ｎは２以上の整数)とし、前記複数の材料の体積分率をそれぞれｖB1，ｖB2，………ｖBnとし、前記仮想複合材料の有効材料定数をＣ*としたとき、
前記非線形方程式は下記式で表され、
下記式中のＡ（ＣB1，Ｃ＊），Ａ（ＣB2，Ｃ＊），………，Ａ（ＣBn，Ｃ＊）は、前記仮想複合材料中の前記球状粒子を囲む周囲の材料定数を、算出すべき前記仮想複合材料の有効材料定数として定めることにより得られ、前記材料定数Ｃ＊に関して非線形な式で定まる比例定数である、請求項１に記載の算出方法。
Ｃ*  ＝  ＣA＋ｖB1・（ＣB1−ＣA）・Ａ（ＣB1，Ｃ＊）
＋ｖB2・（ＣB2−ＣA）・Ａ（ＣB2，Ｃ＊）
………＋ｖBn・（ＣBn−ＣA）・Ａ（ＣBn，Ｃ＊）
【請求項３】
前記仮想複合材料の応力およびひずみを、せん断成分と静水圧成分に分解し、前記仮想複合材料中の、前記静水圧成分に対応する変位をｕi（i＝１、２または３）としたとき、
前記変位ｕiは、前記球状粒子の内部において、前記球状粒子の中心を基準とする位置ｘi（i＝１、２または３）に比例し、
前記変位ｕiは、前記球状粒子の外部において、前記球状粒子の中心を基準とする位置ｘi（i＝１、２または３）に比例する項と、前記球状粒子の中心からの距離の３乗に反比例し、前記位置ｘi（i＝１、２または３）に比例する項を有するように定めることにより、前記非線形方程式が用意される、請求項１または２に記載の算出方法。
【請求項４】
前記仮想複合材料の応力およびひずみを、せん断成分と静水圧成分に分解し、前記仮想複合材料中の、前記せん断成分に対応する変位をｕi（i＝１、２または３）としたとき、
前記変位ｕiを、前記球状粒子の内部および外部において、前記球状粒子の中心からの距離の０乗、距離の２乗、距離の−３乗および距離の５乗のそれぞれに比例する項の加算式で定めることにより、前記非線形方程式が用意される、請求項１または２に記載の算出方法。
【請求項５】
前記非線形方程式は、ニュートン・ラフソン法を用いて、解が収束するように解かれる請求項１〜４のいずれか１項に記載の算出方法。
【請求項６】
力学的材料定数が既知の母相中に、力学的材料定数が既知の複数の種類の材料を含む複合材料に関して、いずれか１種類の材料の体積分率をコンピュータが算出する方法であって、
前記複合材料の有効材料定数を実験結果から特定するステップと、
前記母相中に、前記複数の種類の材料のそれぞれが球状粒子として分散した仮想複合材料を、前記複合材料として定めることにより、この仮想複合材料中のいずれか１種類の材料の体積分率を未知数とする非線形方程式を用意するステップと、
前記用意した非線形方程式を解いて、前記材料の体積分率を求めるステップと、を有し、
前記非線形方程式は、前記仮想複合材料中の前記複数の種類の材料のそれぞれの前記球状粒子を囲む周囲の材料定数を、前記特定された前記複合材料の有効材料定数として定めることにより得られる再帰型非線型方程式である、算出方法。
【請求項７】
前記母相の材料定数をＣAとし、前記複数の種類の材料の材料定数をＣB1，ＣB2，………ＣBn(ｎは２以上の整数)とし、前記仮想複合材料の有効材料定数をＣ*とし、前記複数の材料の体積分率をそれぞれｖB1，ｖB2，………ｖBnとしたとき、
前記非線形方程式は下記式で表され、かつ、Ａ（ＣB1，Ｃ＊），Ａ（ＣB2，Ｃ＊），………，Ａ（ＣBn，Ｃ＊）は、前記仮想複合材料中の前記球状粒子を囲む周囲の材料定数を、算出すべき前記仮想複合材料の有効材料定数として定めることにより得られ、前記材料定数Ｃ＊に関して非線形な式で定まる比例定数である、請求項６に記載の算出方法。
Ｃ*  ＝  ＣA＋ｖB1・（ＣB1−ＣA）・Ａ（ＣB1，Ｃ＊）
＋ｖB2・（ＣB2−ＣA）・Ａ（ＣB2，Ｃ＊）
………＋ｖBn・（ＣBn−ＣA）・Ａ（ＣBn，Ｃ＊）
【請求項８】
前記仮想複合材料の応力およびひずみを、せん断成分と静水圧成分に分解し、前記仮想複合材料中の、前記静水圧成分に対応する変位をｕi（i＝１、２または３）としたとき、
前記変位ｕiは、前記球状粒子の内部において、前記球状粒子の中心を基準とする位置ｘi（i＝１、２または３）に比例し、
前記変位ｕiは、前記球状粒子の外部において、前記球状粒子の中心を基準とする位置ｘi（i＝１、２または３）に比例する項と、前記球状粒子の中心からの距離の３乗に反比例し、前記位置ｘi（i＝１、２または３）に比例する項を有するように定めることにより、前記非線形方程式が用意される、請求項６または７に記載の算出方法。
【請求項９】
前記仮想複合材料の応力およびひずみを、せん断成分と静水圧成分に分解し、前記仮想複合材料中の、前記せん断成分に対応する変位をｕi（i＝１、２または３）としたとき、
前記変位ｕiを、前記球状粒子の内部および外部において、前記球状粒子の中心からの距離の０乗、距離の２乗、距離の−３乗および距離の５乗のそれぞれに比例する項の加算式で定めることにより、前記非線形方程式が用意される、請求項６または７に記載の算出方法。
【請求項１０】
前記非線形方程式は、２分法を用いて解かれる請求項６〜９のいずれか１項に記載の算出方法。
【請求項１１】
力学的材料定数が既知の母相中に、力学的材料定数が既知の複数の種類の材料を含む複合材料に関して、複合材料の力学的有効材料定数を算出する、コンピュータの実行可能なプログラムを記録した記録メディアであって、
前記プログラムは、
前記母相中に、前記複数の種類の材料のそれぞれが既知の体積分率で球状粒子として分散した仮想複合材料を、前記複合材料として定めることにより、この仮想複合材料の有効材料定数を未知数とする非線形方程式を用意する手順と、
前記用意した非線形方程式を解いて、前記仮想複合材料の有効材料定数を、前記複合材料の有効材料定数として求める手順とを、コンピュータに実行させ、
前記非線形方程式は、前記仮想複合材料中の前記複数の種類の材料のそれぞれに対応する前記球状粒子を囲む周囲の材料定数を、算出すべき前記仮想複合材料の有効材料定数として定めることにより得られる再帰型非線型方程式である、記録メディア。
【請求項１２】
力学的材料定数が既知の母相中に、力学的材料定数が既知の複数の種類の材料を含む複合材料に関して、いずれか１種類の材料の体積分率を算出する、コンピュータの実行可能なプログラムを記録した記録メディアであって、
前記プログラムは、
前記複合材料の有効材料定数を実験結果から特定する手順と、
前記母相中に、前記複数の種類の材料のそれぞれが球状粒子として分散した仮想複合材料を、前記複合材料として定めることにより、この仮想複合材料中のいずれか１種類の材料の体積分率を未知数とする非線形方程式を用意する手順と、
前記用意した非線形方程式を解いて、前記材料の体積分率を求める手順とを、コンピュータに実行させ、
前記非線形方程式は、前記仮想複合材料中の前記複数の種類の材料のそれぞれの前記球状粒子を囲む周囲の材料定数を、前記特定された前記複合材料の有効材料定数として定めることにより得られる再帰型非線型方程式である、記録メディア。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、力学的材料定数が既知の母相中に、力学的材料定数が既知の複数種類の材料を含む複合材料に関して、複合材料の力学的有効材料定数を算出する方法と、力学的材料定数が既知の母相中に、力学的材料定数が既知の複数の種類の材料を含む複合材料に関して、いずれか１種類の材料の体積分率を算出する方法と、これらの方法をコンピュータに実行させるプログラムを記録した記録メディアに関する。
【背景技術】
【０００２】
従来より、所定の材料が母相中に分散した複合材料の力学特性を精度よく推定することが盛んに試みられている。これは、複合材料が所望の特性になるように、どのような力学特性を持つ材料を、どの程度の体積分率で分散させればよいか、実験上で探索する替わりに、コンピュータ上で効率よく探索するためである。この探索の結果、所望の材料に関する早期の配合設計が可能となる。
【０００３】
このような状況下、下記特許文献１には、複数の微小要素からなり３次元的に不均質な変形特性を有するミクロ構造が１方向にのみ周期的に配置されるマクロ構造体の構造解析方法が記載されている。当該文献では、マクロ構造体における周期性の単位であるユニットセルを特定し、ユニットセルの均質な材料特性を有するものとして均質化弾性係数を求める。この後、マクロ構造体を、均質化弾性特性を有するものとしてモデル化し、周期性配置の１方向に沿った任意の位置における変形を求めるマクロスケール解析を行う。得られた周期性配置の１方向に沿った任意の位置における変形を、その位置におけるユニットセルを構成する各微小要素のそれぞれに与えてそれらの局所的応答を得るローカル解析を行う。
この構造解析方法により、断面内で不均質なマクロ構造の構造計算時間を短縮することができるとされている。
【０００４】
しかし、上記構造解析方法は、微小要素を用いて構成した有限要素モデルを用いる方法である。このため、上記方法は、モデルの作成及び計算に多くの時間を要し、迅速な初期設計及び初期開発において、有効な手段となりにくい、といった問題があった。
【０００５】
一方、複合材料の力学特性を求めるために、従来より、バネ及びダッシュポットを用いた古典的解析モデルも用いられる。このモデルは、計算時間が短く効率的であるが、複合材料のミクロ状態を考慮したモデルでないため、得られる結果の情報が少なく、精度も低い、といった問題もある。
【先行技術文献】
【特許文献】
【０００６】
【特許文献１】特開２００７−１２２２４２号公報
【発明の概要】
【発明が解決しようとする課題】
【０００７】
そこで、本発明は、上記問題点を解決するために、従来の有限要素モデルを用いた計算方法に比べて効率よく短時間に計算できる複合材料の材料定数の算出方法および複合材料中の材料の体積分率の算出方法を提供するとともに、上記方法をコンピュータで実行するプログラムを記録した記録メディアを提供することを目的とする。
【課題を解決するための手段】
【０００８】
上記目的は、力学的材料定数が既知の母相中に、力学的材料定数が既知の複数の種類の材料を含む複合材料に関して、複合材料の力学的有効材料定数をコンピュータが算出する、以下に示すステップを有する算出方法で達成することができる。
すなわち、この算出方法は、
（Ａ）前記母相中に、前記複数の種類の材料のそれぞれが既知の体積分率で球状粒子として分散した仮想複合材料を、前記複合材料として定めることにより、この仮想複合材料の有効材料定数を未知数とする非線形方程式を用意するステップと、
（Ｂ）前記用意した非線形方程式を解いて、前記仮想複合材料の有効材料定数を、前記複合材料の有効材料定数として求めるステップとを有する。
（Ｃ）その際、前記非線形方程式は、前記仮想複合材料中の前記複数の種類の材料のそれぞれに対応する前記球状粒子を囲む周囲の材料定数を、算出すべき前記仮想複合材料の有効材料定数として定めることにより得られる再帰型非線型方程式である。
【０００９】
また、上記目的は、力学的材料定数が既知の母相中に、力学的材料定数が既知の複数の種類の材料を含む複合材料に関して、いずれか１種類の材料の体積分率をコンピュータが算出する、以下に示すステップを有する算出方法により達成することができる。
すなわち、この算出方法は、
（Ｄ）前記複合材料の有効材料定数を実験結果から特定するステップと、
（Ｅ）前記母相中に、前記複数の種類の材料のそれぞれが球状粒子として分散した仮想複合材料を、前記複合材料として定めることにより、この仮想複合材料中のいずれか１種類の材料の体積分率を未知数とする非線形方程式を用意するステップと、
（Ｆ）前記用意した非線形方程式を解いて、前記材料の体積分率を求めるステップと、を有する。
（Ｇ）その際、前記非線形方程式は、前記仮想複合材料中の前記複数の種類の材料のそれぞれの前記球状粒子を囲む周囲の材料定数を、前記特定された前記複合材料の有効材料定数として定めることにより得られる再帰型非線型方程式である。
【００１０】
さらに、上記目的は、力学的材料定数が既知の母相中に、力学的材料定数が既知の複数の種類の材料を含む複合材料に関して、複合材料の力学的有効材料定数を、上述した算出方法により算出する、コンピュータの実行可能なプログラムを記録した記録メディアによって、達成することができる。
同様に、力学的材料定数が既知の母相中に、力学的材料定数が既知の複数の種類の材料を含む複合材料に関して、いずれか１種類の材料の体積分率を、上述した算出方法により算出する、コンピュータの実行可能なプログラムを記録した記録メディアによって、達成することができる。
【００１１】
上述の算出方法に共通する再帰型非線型方程式は、後述する式（６）で表される方程式に基づく。具体的には、再帰型非線形方程式に用いる式（６）中の比例定数ＡU，ＡFは、仮想複合材料の応力場、ひずみ場をナビエの方程式を解くことにより定まる。そのとき、仮想複合材料中の球状粒子を囲む周囲の材料定数を、仮想複合材料の有効材料定数として定める。すなわち、セルフコンシステント近似を用いる。
【発明の効果】
【００１２】
本発明は、母相中に複数の種類の材料が球状粒子として分散した仮想複合材料を定め、さらに、用いる非線形方程式は、仮想複合材料中のそれぞれの種類の材料に対応する球状粒子を囲む周囲の材料定数を、算出すべき前記仮想複合材料の有効材料定数として定めることにより得られた、解析的な再帰型非線型方程式である。このため、複合材料の有効材料定数の算出および複合材料中の所定の材料の体積分率の算出を効率よく短時間に行うことができる。
【図面の簡単な説明】
【００１３】
【図１】本発明の複合材料の材料定数の算出方法および複合材料中の材料の体積分率の算出方法を実施する計算装置の一例の概略の構成を示す図である。
【図２】（ａ）〜（ｃ）は複合材料の応力ひずみを説明する図である。
【図３】本発明の複合材料の材料定数の算出方法および複合材料中の材料の体積分率の算出方法で用いる均一な材質の仮想複合材料を説明する図である。
【図４】（ａ）及び（ｂ）は、本発明の複合材料の材料定数の算出方法および複合材料中の材料の体積分率の算出方法で用いる仮想複合材料の構造を説明する図である。
【図５】本発明の複合材料の材料定数の算出方法の一例のフローを示す図である。
【図６】本発明の複合材料中の材料の体積分率の算出方法の一例のフローを示す図である。
【図７】複合材料の有限要素モデルを説明する図である。
【発明を実施するための形態】
【００１４】
以下、本発明の複合材料の材料定数の算出方法、および複合材料中の材料の体積分率の算出方法を、下記図面に示す実施形態に基いて詳細に説明する。
図１は、複合材料の材料定数の算出方法、および複合材料中の材料の体積分率の算出方法を実施する計算装置１０を示す。
【００１５】
計算装置１０は、ＣＰＵ１２、ＲＯＭ１４、ＲＡＭ１６および入出力ポート１８を備えるコンピュータにより構成されている。計算装置１０は、ＲＯＭ１４に記録されたプログラムを起動することにより、条件設定モジュール２０と、非線形方程式設定モジュール２２と、方程式解法モジュール２４と、収束解判定モジュール２６と、結果処理モジュール２８と、を形成する。すなわち、上記モジュールは、ソフトウェアを起動することによって構成される各部分である。計算装置１０のモジュールは、例えば、数式処理ソフトウェアを用いて形成される。
計算装置１０は、入出力ポート１８を介して、記憶装置３０、マウスやキーボード等の入力装置３２、および、プリンタやモニタ等の出力装置３４に接続されている。
【００１６】
計算装置１０は、
（１）力学的材料定数が既知の母相中に、力学的材料定数が既知の、複数の種類の材料を含む複合材料に関して、複合材料の力学的有効材料定数を算出する第１の処理と、
（２）力学的材料定数が既知の母相中に、力学的材料定数が既知の、複数の種類の材料を含む複合材料に関して、少なくとも１種類の材料の体積分率を算出する第２の処理と、を選択的に実行する。
ここで、複合材料の力学的有効材料定数とは、複合材料全体の力学材料定数である。複合材料全体の力学材料定数とは、複合材料を均質な材料としてみた時の力学材料定数をいう。以降、力学的材料定数は材料定数という。
【００１７】
条件設定モジュール２０は、第１の処理及び第２の処理のどちらの処理を行うかを設定するとともに、第１の処理あるいは第２の処理に必要な情報を揃えて値を設定する部分である。
このモジュール２０は、第１の処理では、母相の材料定数と、各材料の材料定数、および各材料の体積分率の各数値を設定する。材料定数は、例えば、ラメ定数λ，μの組、あるいは、ラメ定数λ，μの一方と体積弾性率Ｋの組を含む。以降では、ラメ定数λ，μの組を中心に説明するが、ラメ定数λ，μの一方と体積弾性率Ｋの組を用いることも効率的である。
また、このモジュール２０は、第２の処理では、実験により取得した複合材の有効材料定数、各材料の材料定数、および上記材料のうちのいくつかの体積分率の数値を設定する。
これらの数値は、入力装置３２を介したオペレータの指示入力により、あるいは記憶装置３０に予め記憶された情報を呼び出すことにより、設定される。
【００１８】
非線形方程式設定モジュール２２は、条件設定モジュール２０で設定された材料定数および体積分率の値を用いて、第１の処理あるいは第２の処理に応じて、非線形方程式を用意する部分である。このモジュールは、第１の処理及び第２の処理のいずれの場合も、複数の種類の材料が母相中に球状粒子として分散した複合材料を想定したとき、この複合材料の有効材料定数を算出するための再帰型非線形方程式、すなわち、後述する式（６）を呼び出し、この方程式の各係数に、材料定数および体積分率に基づいた値を与える。これにより、解くべき未知数からなる方程式が用意される。
【００１９】
方程式解法モジュール２４は、第１の処理の場合、算出すべき複合材料の有効材料定数に値を付与したとき、この値の付与により、用意された再帰型非線形方程式を用いたニュートン・ラフソン法の漸化式に従って、収束解に近い複合材の材料定数の値を再帰的に算出する部分である。すなわち、方程式解法モジュール２４は、最初に複合材料の有効材料定数の値として初期値（ｎ＝１）を設定し、この値を用いて、漸化式にしたがって、収束解に近い複合材料の有効材料定数の値（ｎ＝２）を再帰的に求める。この値は、後述する収束解判定モジュール２６に送られ、判定の結果、解が収束しない場合、算出された材料定数の値を用いて、上記漸化式に従って、さらに収束解に近い複合材料の有効材料定数の値（ｎ＝３）を再帰的に求める。このようにして、方程式解法モジュール２４は、収束解判定モジュール２６において、収束解が見出されるまで、複合材料の有効材料定数の値を再帰的に求める処理を繰り返す。
【００２０】
また、第２の処理の場合、方程式解法モジュール２４は、下記に示す関数ｆ（ｘ）（ｘは求めようとする所定の材料の体積分率）を定め、この関数ｆ（ｘ）の上限値ｘ＝ｘ1における値ｆ（ｘ1）と下限値ｘ＝ｘ2における値ｆ（ｘ2）との間で積ｆ（ｘ1）・ｆ（ｘ2）を算出し、この結果を収束解判定モジュール２６に送る。
【００２１】
関数ｆ（ｘ）＝（実験により取得した複合材料の有効材料定数）−（再帰型非線形方程式（６）で求める体積分率ｘにおける複合材料の有効材料定数）
【００２２】
このような積ｆ（ｘ1）・ｆ（ｘ2）の算出は、収束解が見出されるまで、上限値ｘ＝ｘ1と下限値ｘ＝ｘ2における値を変更しながら繰り返し求められる。
【００２３】
収束解判定モジュール２６は、方程式解法モジュール２４の算出結果に基づいて収束解か否かを判定する部分である。
第１の処理の場合、繰返回数（ｎ＋１）で求めた複合材料の有効材料定数の値と、繰返回数ｎでもとめた複合材料の有効材料定数の値との差分の絶対値が予め設定された閾値を下回る場合、繰返回数（ｎ＋１）における複合材料の有効材料定数の値が収束解であるとされる。それ以外は、収束解は見出されないとして、繰返回数（ｎ＋１）で求めた値を用いて方程式解法モジュール２４で計算が行われるように、指示される。上述したように、再帰型非線型方程式は、係数等に具体的な数値が付与された方程式なので、この非線形方程式の微分した関数も既知である。したがって、収束解判定モジュール２６は、公知のニュートン・ラフソン法を用いて収束解を見出すことができる。
第２の処理の場合、ｆ（ｘ1）・ｆ（ｘ2）の値が、負か否かが調べられ、後述するように、２分法を用いて、上限値および下限値の値が再設定される。
【００２４】
結果処理モジュール２８は、収束解判定モジュール２６で得られた収束解を用いて、複合材料の有効ヤング率やの有効せん断剛性を算出し、あるいは、所定の材料の体積分率を取得する部分である。
【００２５】
出力装置３４は、結果処理モジュール２８で得られた各情報をプリント出力し、あるいは画面表示を行う部分である。
記憶装置３０には、予め設定された複合材料や複合材料に含まれる材料の材料定数等が記録保持されているデータベースを備えている。
【００２６】
なお、上述の方程式解法モジュール２４は、母相αに２つの異なる材料β，γを含む場合、母相αに材料β、γが分散した均一の材質の仮想複合材料を定めて、この仮想複合材料の有効材料定数を算出する。
【００２７】
このような処理は、複合材料の有効材料定数の算出を、以下に示す再帰型非線形方程式を用いて解析的に行うことにより実現される。以下、この再帰型非線形方程式について詳細に説明する。
【００２８】
図２（ａ）〜（ｃ）は、複合材料を説明する図である。
複合材料５０は、図２（ａ）に示すように、エポキシ樹脂（以降、エポキシという）Ｐからなる母相中に、ウレタンＵ及びフィラーＦを含む。
図２（ｂ）および（ｃ）は、複合材料５０に作用する応力およびひずみを示している。  以下で説明する応力、ひずみおよび材料定数は、それぞれ２階テンソル、２階テンソル、および４階テンソルであるが、わかり易く説明するためにスカラーの形式で説明する。
複合材料５０では、図２（ｂ）に示すように、ウレタンＵに作用する応力σUと、フィラーＦに作用する応力σFと、エポキシＰに作用する応力σPとが加算されて、複合材料全体の平均応力σTとなる。これに対応して、図２（ｃ）に示すように、ウレタンＵに作用するひずみεUと、フィラーＦに作用するひずみεFと、エポキシＰに作用するひずみεPとが加算されて、複合材料全体の平均ひずみεTとなる。平均応力σTおよび平均ひずみεTは、下記式（１）、（２）で表される。平均応力及び平均ひずみは、いずれも体積分率を重み付け係数とした体積平均の応力及び体積平均のひずみを意味する。
【００２９】
【数１】
2010205254.tif 000003
【数２】
2010205254.tif 000004
【００３０】
一方、複合材料全体における平均応力σTとひずみεTとの関係を規定する材料定数Ｃ＊は、図３（ａ）に示すように、下記式（３）で表される。すなわち、材料定数Ｃ＊を持つ均一な材質の仮想複合材料が想定される。
【数３】
2010205254.tif 000005
【００３１】
ここで、ひずみに注目し、ウレタンＵに作用するひずみεUと複合材料５０全体に作用するひずみεTとの間の関係を、下記式（４）で表す。同様に、フィラーＦに作用するひずみεFと複合材料５０全体に作用するひずみεTとの間の関係を、下記式（５）で表す。
【数４】
2010205254.tif 000006
【数５】
2010205254.tif 000007
【００３２】
式（４）は、ウレタンＵの平均ひずみであるεUがεTの関数として表され、εUが比例定数ＡUによりεTと関係付けられていることを表している。同様に、式（５）は、フィラーＦの平均ひずみであるεFがεTの関数として表され、εFが比例定数ＡFによりεTと関係付けられていることを表している。
この比例定数ＡUおよびＡFは、後述するように、ウレタンＵおよびフィラーＦを球状粒子と考え、この球状粒子を囲む周囲が、均一な材質の仮想複合材料であると考えることにより求めることができる。すなわち、比例定数ＡUを、ウレタンＵの材料定数と均一な材質の仮想複合材料の有効材料定数Ｃ＊とを用いて表すことができる。比例定数ＡFを、フィラーＦの材料定数と均一な材質の仮想複合材料の有効材料定数Ｃ＊とを用いて表すことができる。ウレタンＵの材料定数をＣUとすると、ＡU＝Ａ（ＣU，Ｃ＊）と表される。フィラーＦの材料定数をＣFとすると、ＡF＝Ａ（ＣF，Ｃ＊）と表される。ＡU＝Ａ（ＣU，Ｃ＊）およびＡF＝Ａ（ＣF，Ｃ＊）の求め方は後述する。Ａ（ＣU，Ｃ＊）およびＡ（ＣF，Ｃ＊）はいずれも、ＣU，Ｃ＊およびＣF，Ｃ＊の複雑な式で表された非線形な式となっている。ここで、ＣU，ＣF，Ｃ＊は、例えば、それぞれの材料のラメ定数λU、μU、λF、μFおよびλ＊，μ＊  を表す。このような比例定数ＡU，ＡFを用いて、式（１）〜（５）を整理すると、下記式（６）が導かれる。
【００３３】
【数６】
2010205254.tif 000008
【００３４】
式（６）は、左辺の材料定数Ｃ＊を求める式であり、右辺にある比例定数ＡU，ＡFはそれぞれ、ＣU，Ｃ＊およびＣF，Ｃ＊に関して非線形な式となっている。このため、式（６）は、材料定数Ｃ＊に関して、再帰型の非線形方程式となっている。
【００３５】
このような比例定数ＡU，ＡFは、図４（ａ），（ｂ）に示すように、材料定数Ｃ＊を持つ均一な材質の仮想複合材料中に球状粒子であるウレタンＵおよびフィラーＦが配置されるモデルを仮想し、この球状粒子を囲む周囲の材料定数を、算出すべき仮想複合材料の有効材料定数Ｃ＊として定めることにより得られる。
すなわち、比例定数ＡU，ＡFは、体積力が作用せず、無限遠で印加されるひずみε∞ij（ｉ，ｊは、１〜３の自然数）が存在する条件下、公知のナビエの方程式を解くことにより得られる。なお、無限遠で印加されるひずみε∞ijは下記式（７）のように、第１項の静水圧成分と第２項のせん断成分とに分解できる。このため、比例定数ＡU，ＡFについても静水圧成分と第２項のせん断成分に分けて算出することができる。比例定数ＡU，ＡFは、いずれも球状粒子の材料定数がＣUかＣFの違いしかないので、以下の説明では、球状粒子の材料定数をＣU，ＣFの代表としてＣBを用いて説明する。
【００３６】
【数７】
2010205254.tif 000009
【００３７】
以下、Ａ（ＣB，Ｃ＊）の第１項の静水圧成分と第２項のせん断成分の算出について説明する。
【００３８】
（静水圧成分に基づくＡ（ＣB，Ｃ＊）の算出）
式（７）中の第１項の静水圧成分に関するナビエの方程式を満足する変位ｕi（ｉは１〜３の自然数）は、下記式（８）のように表される。これは、右辺のテンソルの次元を変位ｕiの次元と対応させる必要性から一意的に定まる。
ここで、球状粒子における変位ｕiが原点で有限の値を持ち、ひずみが無限遠で式（７）中の第１項のひずみになる条件と、変位ｕiおよび球体粒子の表面における表面張力の連続性の条件とを与えることにより、球体粒子内部の変位ｕiと球体粒子外部の変位ｕiを求めることができる。具体的には、球体粒子内部の変位ｕiは下記式（９）で、球体粒子外部の変位ｕiは下記式（１０）で表される。
【００３９】
【数８】
2010205254.tif 000010
【数９】
2010205254.tif 000011
【数１０】
2010205254.tif 000012
【００４０】
すなわち、母相に球状粒子がある空間内の各方向の変位ｕiは、球状粒子の内部において、球状粒子の中心を基準として定められる位置ｘi（i＝１、２または３）に比例する。変位ｕiは、さらに、球状粒子の外部において、球状粒子の中心を基準として定められる位置ｘi（i＝１、２または３）に比例する項と、球状粒子の中心からの距離の３乗に反比例し、位置ｘiに比例する項とを有するように定められる。これにより、比例定数Ａ（ＣB，Ｃ＊）の静水圧成分を求めることができる。
【００４１】
式（９）及び式（１０）から、球状粒子内部における応力σijは下記式（１１）で表すことができ、球状粒子内部におけるひずみεijは、下記式（１２）で表すことができる。この式（１２）は、上述した式（５）の関係式の静水圧成分に対応するものである。したがって、Ａ（ＣB，Ｃ＊）の静水圧成分は、式（１３）のように表される。したがって、ウレタンＵの比例定数ＡUは、式（１３）中のλBおよびμBをλUおよびμUに置き換えるとよい。同様に、フィラーＦの比例定数ＡFは、式（１３）中のλBおよびμBをλFおよびμFに置き換えるとよい。
【００４２】
【数１１】
2010205254.tif 000013
【数１２】
2010205254.tif 000014
【数１３】
2010205254.tif 000015
【００４３】
（せん断成分に基づくＡ（ＣB，Ｃ＊）の算出）
式（７）中の第２項のせん断成分に関するナビエの方程式を満足する変位ｕi（ｉは１〜３の自然数）は、球状粒子の内部及び外部を問わず、球状粒子の中心からの距離の０乗、距離の２乗、距離の−３乗および距離の５乗のそれぞれに比例する項の加算式で定められる。具体的には、変位ｕiを、｛ｆ（ｒ）・（ｘiｘj／ｒ2）＋ｇ（ｒ）・δij｝・ｘk（ｉ，ｊ，ｋ＝１〜３の整数）に比例するように定める。この式は、テンソルの階数をナビエの方程式上で対応させる必要性から一意的に定められたものである。ここで、ｘｊ，ｘkは、３次元座標におけるｘiと異なる座標成分を表す。また、ｆ（ｒ）、ｇ（ｒ）は球状粒子の中心からの距離ｒのみによる多項式を表す。このとき、ナビエの方程式を満足するｆ（ｒ）、ｇ（ｒ）をｆ（ｒ）∝ｒｌ、ｇ（ｒ）∝ｒｌ（ｌは整数）としてｒの次数ｌを次元解析により求めることにより、ｒの次数ｌは０，２，−３，５となる。したがって、ナビエの方程式を満足する変位ｕiは、球状粒子の内部及び外部を問わず、球状粒子の中心からの距離の０乗、距離の２乗、距離の−３乗および距離の５乗のそれぞれに比例する項の加算式により表される。例えば、次数ｌ＝２のときの式は、下記式（１４）に示すように表される。勿論、式（１４）中のλ、μは、球状粒子の内部では、球状粒子の内部の材料定数が、球状粒子の外部では、球状粒子の外部の材料定数が用いられる。
【数１４】
2010205254.tif 000016
【００４４】
すなわち、変位ｕiは、球状粒子の内部および外部において、球状粒子の中心からの距離の０乗、距離の２乗、距離の−３乗および距離の５乗のそれぞれに比例する項の加算式として表すことができる。これにより、再帰型非線形方程式を得るための比例定数Ａ（ＣB，Ｃ＊）のせん断成分を求めることができる。
【００４５】
具体的には、球状粒子の表面における変位ｕiおよび球体粒子の表面における表面張力の連続性の条件を与えて、下記式（１５）に示す球体粒子内部の応力が求められる。この式（１５）から、比例定数Ａ（ＣB，Ｃ＊）のせん断成分が下記式（１６）のように表される。したがって、ウレタンＵの比例定数ＡUは、式（１６）中のμBをμUに置き換えるとよい。同様に、フィラーＦの比例定数ＡFは、式（１６）中のμBをμFに置き換えるとよい。
【数１５】
2010205254.tif 000017
【数１６】
2010205254.tif 000018
【００４６】
以上より、ウレタンＵの比例定数ＡUおよびフィラーＦの比例定数ＡFが定められ、これらの比例定数が、材料定数Ｃ＊を求める式（６）に代入されて、再帰型非線形方程式が導かれる。この再帰型非線形方程式は、ウレタンＵおよびフィラーＦが球状粒子としてエポキシＰに分散された均質な材質の仮想複合材料を仮想し、この仮想複合材料中の球状粒子を囲む周囲の材料定数を、算出すべき仮想複合材料の有効材料定数として定めることにより実現したものである。
この仮想複合材料を利用した本発明の算出方法による複合材料の有効材料定数の計算結果は、後述するように、従来の長時間の計算を要する有限要素モデルを用いた場合に近く、精度が高いものとなっている。
【００４７】
このような再帰型非線形方程式を用いて複合材料の有効材料定数Ｃ＊を算出する方法について説明する。
【００４８】
図５は、上記材料定数Ｃ＊を算出するフローを説明する図である。
まず、条件設定モジュール２０にて、複合材料の有効材料定数Ｃ＊を算出する第１の処理が設定される。この設定は、入力装置３２を介してオペレータから指示入力されて行われる。ここで、材料定数Ｃ＊は、算出すべき未知数であるラメの定数λef,μefを代表して表している。
次に、ウレタンＵ、フィラーＦおよびエポキシＰの材料定数ＣU，ＣF，ＣP、すなわち、λP，μP，λU，μU，λF，μFと、ウレタンＵおよびフィラーＦの体積分率ｖu，ｖFの値が、記憶装置３０中のデータベースから呼び出されて取得される（ステップＳ１０）。材料定数と体積分率の値は、非線形方程式設定モジュール２２に送られて、上述の式（６）の各係数の値が定められ、材料定数Ｃ＊のみの再帰型非線形方程式に修正される。すなわち、方程式が用意される。
【００４９】
次に、方程式解法モジュール２４にて、材料定数Ｃ＊の未知数であるラメの定数λef,μefの初期値が設定される（ステップＳ２０）。初期値を設定する理由は、式（６）が非線形であり、解が収束するまで繰り返すためである。初期値はどのような値であってもよいが、例えば、母相であるエポキシＰの材料定数ＣPとフィラーＦの材料定数ＣFの、体積分率による加重平均値を選択するとよい。この初期値を、ラメの定数λef,＝λ1，μef＝μ1と定める。さらに、この初期値をそれぞれ、λn，μn（ｎ＝１）と表す。
次に、方程式解法モジュール２４にて、上述の式（６）を修正した下記式（１７）を定義する。具体的には、式（６）の左辺を右辺に移動し、その時の右辺をＦ（Ｃ＊)として定義する。ＡＵはＡ（ＣＵ，Ｃ＊）であり、ＡＦはＡ（ＣＦ，Ｃ＊）である。
【００５０】
【数１７】
2010205254.tif 000019
【００５１】
さらに、下記式（１８）にしたがって、Ｆ（Ｃ＊)を用いてＣ(n+1)が定義される。ここで、式（１８）は、Ｆ（Ｃ＊)＝０となる収束解を求めるためにニュートン・ラフソン法を適用した式である。ここで、Ｃ(n)，Ｃ(n+1)はラメの定数λn，μn，λn+1，μn+1を代表して表す符号である。式（１８）のＣ(n)に、ラメの定数λn，μnが与えられ、式（１８）に従ってラメの定数λn+1，μn+1が算出され取得される（ステップＳ３０）。ここで、式（１８）中のＦ´（Ｃ(n)）はヤコビのテンソルであり、Ｆ´（Ｃ(n)）-1はＦ´（Ｃ(n)）の逆テンソルを表す。すなわち、Ｆ´（Ｃ(n)）-1・Ｆ’（Ｃ(n)）＝Ｉ（Ｉは単位テンソル）である。
【００５２】
【数１８】
2010205254.tif 000020
【００５３】
算出されたλn+1，μn+1と、この算出に用いられたλn，μnが収束解判定モジュール２６に送られて、それぞれの差分λn+1−λn，μn+1−μnの絶対値が求められる。この絶対値が予め設定された閾値ε1，ε2と比較される（ステップＳ４０）。比較の結果、差分λn+1−λn，μn+1−μnの絶対値がいずれも閾値ε1，ε2より小さい条件を満たす場合（Ｙｅｓの場合）、算出すべき仮想複合材料の有効材料定数Ｃ＊であるラメの定数λ*,μ＊は、それぞれλn+1，μn+1に決定される（ステップＳ５０）。一方、比較の結果、上記条件を満たさない場合（Ｎｏの場合）、ｎ＋１をｎとし（ステップＳ６０）、さらに、式（１８）を用いてλn+1，μn+1を算出するためにステップＳ３０に戻る。こうして、ステップＳ４０における条件を満足するまで、ステップＳ３０，Ｓ４０，Ｓ６０が繰り返される。
【００５４】
次に、このように定められた複合材料の有効材料定数Ｃ＊、具体的にはラメの定数λ*,μ＊を用いてヤング率Ｅ*が算出される。ヤング率Ｅ*は、下記式（１９）に従って算出され、プリンタやディスプレイ等の出力装置に出力される。
【数１９】
2010205254.tif 000021
【００５５】
算出されたヤング率Ｅ*は、後述するように、有限要素モデルを用いて数値計算をおこなって得られるヤング率Ｅとよく対応した値となっている。このため、本発明の複合材料の材料定数の算出方法は有効であることがわかる。また、本発明の複合材料の材料定数の算出方法は、式（６）に示される式を用いて解析的に算出できるので、有限要素モデルを用いた計算に比べて短時間に算出結果を得ることができ、効率的である。
【００５６】
本発明の複合材料の材料定数の算出方法で用いる式（６）は、解析式であるので、例えば、母相に含まれるフィラーＦ等の材料の体積分率を解くべき未知数とすることもできる。このとき、複合材料の有効材料定数は、予め実験等により求めておく必要がある。このような所定の材料の体積分率の算出は、上述した第２の処理に対応する。この第２の処理の具体的な方法を以下説明する。
【００５７】
図６は、第２の処理のフローを示す図である。エポキシＰに、フィラーＦとウレタンＵを含む複合材料において、フィラーＦの体積分率ｖFを算出する場合を想定する。
【００５８】
まず、条件設定モジュール２０において、エポキシＰ、ウレタンＵ、フィラーＦの材料定数、具体的には、λP，μP，λU，μU，λF，μFと、ウレタンＵの体積分率ｖUが取得される（ステップＳ１００）。これらの材料定数等は、記憶装置３０中のデータベースから呼び出されて取得される。あるいは、入力装置３２から入力される。
次に、複合材料の力学的実験による計測結果より、複合材料全体の材料定数である有効材料、具体的にはλ＊，μ＊が特定され取得される（ステップＳ１１０）。計測結果は、予め記憶装置３０に記憶されており、条件設定モジュール２０に計測結果が呼び出されて、複合材料の有効材料定数が特定される。非線形方程式設定モジュール２２では、ステップＳ１００、Ｓ１１０にて取得された材料定数の各値が式（６）の各係数に代入されて、未知数が体積分率ｖFの方程式が用意される。さらに、用意された式（６）を用いて、関数ｆ（ｘ）＝（実験により取得した複合材料の有効材料定数）−（式（６）で求める体積分率ｘにおける複合材料の有効材料定数）が定義される。
【００５９】
この後、フィラーＦの体積分率ｖFの上限値ｘ1及び下限値ｘ2が設定される（ステップＳ１２０）。上限値および下限値は、後述する体積分率ｖFの算出を、２分法を用いて行うために設定される。この設定は、入力装置３２を通して、オペレータの指示入力に基づいて設定されてもよいし、予め設定されたデフォルト値を用いてもよい。
【００６０】
次に、方程式解法モジュール２４において、上限値ｘ1及び下限値ｘ2が用いられて、ｆ（ｘ1）・ｆ（ｘ2）の値が負か否かが判定される（ステップＳ１３０）。ｆ（ｘ1）・ｆ（ｘ2）の値が正の場合、上限値ｘ1及び下限値ｘ2が変更される（ステップＳ１４０）。変更の方法は、特に制限されないが、ｘ1及びｘ2は体積分率であるので、上限値ｘ1はより大きく、下限値ｘ2はより小さくなるように設定されるとよい。
【００６１】
ステップＳ２３０における判定で肯定される場合（Ｙｅｓの場合）、公知の２分法に基づいて、以下の処理が行われる。
すなわち、まず、ｘ3＝（ｘ1＋ｘ2）／２が求められ（ステップＳ１５０）、さらに、ｆ（ｘ1）・ｆ（ｘ3）の値が負か否かが判定される（ステップＳ１６０）。判定の結果、ｆ（ｘ1）・ｆ（ｘ3）が負である場合（Ｙｅｓの場合）、下限値ｘ2がｘ3に変更され（ステップＳ１７０）、ｆ（ｘ1）・ｆ（ｘ3）が負でない場合（Ｎｏの場合）、上限値ｘ1がｘ3に変更される（ステップＳ１８０）。
この後、上限値ｘ1と下限値ｘ2との差分の絶対値が予め設定された閾値ε5より小さいか否かが判定される（ステップＳ１９０）。絶対値が閾値ε5より小さい場合（Ｙｅｓの場合）、上限値ｘ1または下限値ｘ2が算出すべき体積分率ｖFとして定められる（ステップＳ２００）。ステップＳ１９０において、絶対値が閾値ε5より小さくない場合（Ｎｏの場合）、ステップＳ１５０に戻る。こうして、ステップＳ１９０の判定で、絶対値が閾値ε5より小さくなるまで、ステップＳ１５０、Ｓ１６０、Ｓ１７０、Ｓ１８０を繰り返す。
【００６２】
以上の方法により、フィラーＦの体積分率ｖFは算出される。算出したフィラーＦの体積分率ｖFは、出力装置３４に出力される。
なお、算出すべき体積分率は、フィラーＦを対象とする他、ウレタンＵを対象とすることもできる。また、算出すべき体積分率は、複数の種類の材料を対象としてもよい。
複合材料の材料定数の算出方法および体積分率の算出方法では、母相に２つの材料が分散した複合材料を用いて説明したが、３つ以上であってもよい。母相は、エポキシに限られず、金属やセラミック等の無機材料であってもよい。また、母相に含まれる材料は、母相を強化する強化材料に限られず、上述したようにウレタン等の軟質材料であってもよい。
【００６３】
また、複合材料の材料定数の算出方法および体積分率の算出方法は、コンピュータ上でプログラムを実行することにより実現することもできる。このプログラムはコンピュータで読み取り可能な記録メディアに記録されたものである。この記録メディアは、通信回線を利用してダウンロードされたプログラムを記録したものも含まれる。
【００６４】
次に、上述の複合材料の有効材料定数の算出方法の有効性を説明する。
複合材料は、以下に示すものが用いられた。
（１）エポキシ、ウレタン：６０ｇ
ウレタン変性エポキシ：          ３０ｇ
ＮＢＲ変性エポキシ  ：          ２０ｇ
ビスフェノールＡ型液状エポキシ：１０ｇ
（２）フィラー：３５ｇ
シリカ                        ：１５ｇ
炭酸カルシウム                ：１０ｇ
酸化カルシウム                ：１０ｇ
（３）硬化系：５ｇ
ＤＩＣＹ（ジシアンジアミド）  ：５ｇ
【００６５】
複合材料におけるヤング率及びポアソン比と、体積分率（％）は以下の通りである。なお、ヤング率およびポアソン比から、公知の式により、ラメの定数に変換して用いた。
【表１】
2010205254.tif 000022
【００６６】
本発明の複合材料の材料定数の算出方法を用いた複合材料の有効材料定数として得られるヤング率は２．１２６（ＧＰａ）であった。一方、有限要素モデルを用いた複合材料の解析により得られたヤング率は２．００７（ＧＰａ）であった。
これより、本発明の方法により得られるヤング率は、有限要素モデルを用いて得られるヤング率および実験により得られたヤング率と極めて近い値を示し、本発明の方法が有効であることがわかる。
【００６７】
なお、有限要素モデルを用いた材料定数の算出は、図７に示すような、エポキシＰを母相とし、このエポキシＰ中に球状のウレタンＵが体積分率３２％含まれ、フィラーＦが体積分率１４％含まれて構成される複合材料を想定した有限要素モデルを作成して、複合材料の有効材料定数を算出した。図中上端部に引っ張り変位を与えて応力ーひずみの初期傾きを調べてヤング率とした。なお、作成した有限要素モデルは、要素数６５５３６個、節点数６６０４９個の２次元平面応力要素モデルである。解析は、非線形有限要素法ソフトウェアＡＢＡＱＵＳ（商品名）を用いた。
【００６８】
以上、本発明の複合材料の材料定数の算出方法、複合材料中の材料の体積分率の算出方法および記録メディアについて詳細に説明したが、本発明は上記実施形態に限定されず、本発明の主旨を逸脱しない範囲において、種々の改良や変更をしてもよいのはもちろんである。
【符号の説明】
【００６９】
１０  計算装置
１２  ＣＰＵ
１４  ＲＯＭ
１６  ＲＡＭ
１８  入出力ポート
２０  条件設定モジュール
２２  非線形方程式設定モジュール
２４  方程式解法モジュール
２６  収束解判定モジュール
２８  結果処理モジュール
３０  記憶装置
３２  入力装置
【図１】
2010205254.tif 000023
【図２】
2010205254.tif 000024
【図３】
2010205254.tif 000025
【図４】
2010205254.tif 000026
【図５】
2010205254.tif 000027
【図６】
2010205254.tif 000028
【図７】
2010205254.tif 000029