JP2010205254A - 複合材料の力学的材料定数の算出方法、複合材料中の材料の体積分率の算出方法および記録メディア - Google Patents

Abstract

【課題】複合材料の材料定数の算出を、効率よく短時間に算出する。また、複合材料中の材料の体積分率の算出を効率よく短時間に行う。
【解決手段】複合材料について、母相中に、複数の種類の材料が既知の体積分率で球状粒子として分散した仮想複合材料を定めることにより、この仮想複合材料の有効材料定数を未知数とする非線形方程式を用意する。次に、この方程式を解いて、仮想複合材料の有効材料定数を、複合材料の有効材料定数として求める。このとき、方程式は、セルフコンシステント法を用いて得られる再帰型の式である。この方程式を用いて複合材料中の材料の体積分率を算出する。
【選択図】 図１

Description

本発明は、力学的材料定数が既知の母相中に、力学的材料定数が既知の複数種類の材料を含む複合材料に関して、複合材料の力学的有効材料定数を算出する方法と、力学的材料定数が既知の母相中に、力学的材料定数が既知の複数の種類の材料を含む複合材料に関して、いずれか１種類の材料の体積分率を算出する方法と、これらの方法をコンピュータに実行させるプログラムを記録した記録メディアに関する。

従来より、所定の材料が母相中に分散した複合材料の力学特性を精度よく推定することが盛んに試みられている。これは、複合材料が所望の特性になるように、どのような力学特性を持つ材料を、どの程度の体積分率で分散させればよいか、実験上で探索する替わりに、コンピュータ上で効率よく探索するためである。この探索の結果、所望の材料に関する早期の配合設計が可能となる。

このような状況下、下記特許文献１には、複数の微小要素からなり３次元的に不均質な変形特性を有するミクロ構造が１方向にのみ周期的に配置されるマクロ構造体の構造解析方法が記載されている。当該文献では、マクロ構造体における周期性の単位であるユニットセルを特定し、ユニットセルの均質な材料特性を有するものとして均質化弾性係数を求める。この後、マクロ構造体を、均質化弾性特性を有するものとしてモデル化し、周期性配置の１方向に沿った任意の位置における変形を求めるマクロスケール解析を行う。得られた周期性配置の１方向に沿った任意の位置における変形を、その位置におけるユニットセルを構成する各微小要素のそれぞれに与えてそれらの局所的応答を得るローカル解析を行う。
この構造解析方法により、断面内で不均質なマクロ構造の構造計算時間を短縮することができるとされている。

しかし、上記構造解析方法は、微小要素を用いて構成した有限要素モデルを用いる方法である。このため、上記方法は、モデルの作成及び計算に多くの時間を要し、迅速な初期設計及び初期開発において、有効な手段となりにくい、といった問題があった。

一方、複合材料の力学特性を求めるために、従来より、バネ及びダッシュポットを用いた古典的解析モデルも用いられる。このモデルは、計算時間が短く効率的であるが、複合材料のミクロ状態を考慮したモデルでないため、得られる結果の情報が少なく、精度も低い、といった問題もある。

特開２００７−１２２２４２号公報

そこで、本発明は、上記問題点を解決するために、従来の有限要素モデルを用いた計算方法に比べて効率よく短時間に計算できる複合材料の材料定数の算出方法および複合材料中の材料の体積分率の算出方法を提供するとともに、上記方法をコンピュータで実行するプログラムを記録した記録メディアを提供することを目的とする。

上記目的は、力学的材料定数が既知の母相中に、力学的材料定数が既知の複数の種類の材料を含む複合材料に関して、複合材料の力学的有効材料定数をコンピュータが算出する、以下に示すステップを有する算出方法で達成することができる。
すなわち、この算出方法は、
（Ａ）前記母相中に、前記複数の種類の材料のそれぞれが既知の体積分率で球状粒子として分散した仮想複合材料を、前記複合材料として定めることにより、この仮想複合材料の有効材料定数を未知数とする非線形方程式を用意するステップと、
（Ｂ）前記用意した非線形方程式を解いて、前記仮想複合材料の有効材料定数を、前記複合材料の有効材料定数として求めるステップとを有する。
（Ｃ）その際、前記非線形方程式は、前記仮想複合材料中の前記複数の種類の材料のそれぞれに対応する前記球状粒子を囲む周囲の材料定数を、算出すべき前記仮想複合材料の有効材料定数として定めることにより得られる再帰型非線型方程式である。

また、上記目的は、力学的材料定数が既知の母相中に、力学的材料定数が既知の複数の種類の材料を含む複合材料に関して、いずれか１種類の材料の体積分率をコンピュータが算出する、以下に示すステップを有する算出方法により達成することができる。
すなわち、この算出方法は、
（Ｄ）前記複合材料の有効材料定数を実験結果から特定するステップと、
（Ｅ）前記母相中に、前記複数の種類の材料のそれぞれが球状粒子として分散した仮想複合材料を、前記複合材料として定めることにより、この仮想複合材料中のいずれか１種類の材料の体積分率を未知数とする非線形方程式を用意するステップと、
（Ｆ）前記用意した非線形方程式を解いて、前記材料の体積分率を求めるステップと、を有する。
（Ｇ）その際、前記非線形方程式は、前記仮想複合材料中の前記複数の種類の材料のそれぞれの前記球状粒子を囲む周囲の材料定数を、前記特定された前記複合材料の有効材料定数として定めることにより得られる再帰型非線型方程式である。

さらに、上記目的は、力学的材料定数が既知の母相中に、力学的材料定数が既知の複数の種類の材料を含む複合材料に関して、複合材料の力学的有効材料定数を、上述した算出方法により算出する、コンピュータの実行可能なプログラムを記録した記録メディアによって、達成することができる。
同様に、力学的材料定数が既知の母相中に、力学的材料定数が既知の複数の種類の材料を含む複合材料に関して、いずれか１種類の材料の体積分率を、上述した算出方法により算出する、コンピュータの実行可能なプログラムを記録した記録メディアによって、達成することができる。

上述の算出方法に共通する再帰型非線型方程式は、後述する式（６）で表される方程式に基づく。具体的には、再帰型非線形方程式に用いる式（６）中の比例定数Ａ^U，Ａ^Fは、仮想複合材料の応力場、ひずみ場をナビエの方程式を解くことにより定まる。そのとき、仮想複合材料中の球状粒子を囲む周囲の材料定数を、仮想複合材料の有効材料定数として定める。すなわち、セルフコンシステント近似を用いる。

本発明は、母相中に複数の種類の材料が球状粒子として分散した仮想複合材料を定め、さらに、用いる非線形方程式は、仮想複合材料中のそれぞれの種類の材料に対応する球状粒子を囲む周囲の材料定数を、算出すべき前記仮想複合材料の有効材料定数として定めることにより得られた、解析的な再帰型非線型方程式である。このため、複合材料の有効材料定数の算出および複合材料中の所定の材料の体積分率の算出を効率よく短時間に行うことができる。

本発明の複合材料の材料定数の算出方法および複合材料中の材料の体積分率の算出方法を実施する計算装置の一例の概略の構成を示す図である。 （ａ）〜（ｃ）は複合材料の応力ひずみを説明する図である。 本発明の複合材料の材料定数の算出方法および複合材料中の材料の体積分率の算出方法で用いる均一な材質の仮想複合材料を説明する図である。 （ａ）及び（ｂ）は、本発明の複合材料の材料定数の算出方法および複合材料中の材料の体積分率の算出方法で用いる仮想複合材料の構造を説明する図である。 本発明の複合材料の材料定数の算出方法の一例のフローを示す図である。 本発明の複合材料中の材料の体積分率の算出方法の一例のフローを示す図である。 複合材料の有限要素モデルを説明する図である。

以下、本発明の複合材料の材料定数の算出方法、および複合材料中の材料の体積分率の算出方法を、下記図面に示す実施形態に基いて詳細に説明する。
図１は、複合材料の材料定数の算出方法、および複合材料中の材料の体積分率の算出方法を実施する計算装置１０を示す。

計算装置１０は、ＣＰＵ１２、ＲＯＭ１４、ＲＡＭ１６および入出力ポート１８を備えるコンピュータにより構成されている。計算装置１０は、ＲＯＭ１４に記録されたプログラムを起動することにより、条件設定モジュール２０と、非線形方程式設定モジュール２２と、方程式解法モジュール２４と、収束解判定モジュール２６と、結果処理モジュール２８と、を形成する。すなわち、上記モジュールは、ソフトウェアを起動することによって構成される各部分である。計算装置１０のモジュールは、例えば、数式処理ソフトウェアを用いて形成される。
計算装置１０は、入出力ポート１８を介して、記憶装置３０、マウスやキーボード等の入力装置３２、および、プリンタやモニタ等の出力装置３４に接続されている。

計算装置１０は、
（１）力学的材料定数が既知の母相中に、力学的材料定数が既知の、複数の種類の材料を含む複合材料に関して、複合材料の力学的有効材料定数を算出する第１の処理と、
（２）力学的材料定数が既知の母相中に、力学的材料定数が既知の、複数の種類の材料を含む複合材料に関して、少なくとも１種類の材料の体積分率を算出する第２の処理と、を選択的に実行する。
ここで、複合材料の力学的有効材料定数とは、複合材料全体の力学材料定数である。複合材料全体の力学材料定数とは、複合材料を均質な材料としてみた時の力学材料定数をいう。以降、力学的材料定数は材料定数という。

条件設定モジュール２０は、第１の処理及び第２の処理のどちらの処理を行うかを設定するとともに、第１の処理あるいは第２の処理に必要な情報を揃えて値を設定する部分である。
このモジュール２０は、第１の処理では、母相の材料定数と、各材料の材料定数、および各材料の体積分率の各数値を設定する。材料定数は、例えば、ラメ定数λ，μの組、あるいは、ラメ定数λ，μの一方と体積弾性率Ｋの組を含む。以降では、ラメ定数λ，μの組を中心に説明するが、ラメ定数λ，μの一方と体積弾性率Ｋの組を用いることも効率的である。
また、このモジュール２０は、第２の処理では、実験により取得した複合材の有効材料定数、各材料の材料定数、および上記材料のうちのいくつかの体積分率の数値を設定する。
これらの数値は、入力装置３２を介したオペレータの指示入力により、あるいは記憶装置３０に予め記憶された情報を呼び出すことにより、設定される。

非線形方程式設定モジュール２２は、条件設定モジュール２０で設定された材料定数および体積分率の値を用いて、第１の処理あるいは第２の処理に応じて、非線形方程式を用意する部分である。このモジュールは、第１の処理及び第２の処理のいずれの場合も、複数の種類の材料が母相中に球状粒子として分散した複合材料を想定したとき、この複合材料の有効材料定数を算出するための再帰型非線形方程式、すなわち、後述する式（６）を呼び出し、この方程式の各係数に、材料定数および体積分率に基づいた値を与える。これにより、解くべき未知数からなる方程式が用意される。

方程式解法モジュール２４は、第１の処理の場合、算出すべき複合材料の有効材料定数に値を付与したとき、この値の付与により、用意された再帰型非線形方程式を用いたニュートン・ラフソン法の漸化式に従って、収束解に近い複合材の材料定数の値を再帰的に算出する部分である。すなわち、方程式解法モジュール２４は、最初に複合材料の有効材料定数の値として初期値（ｎ＝１）を設定し、この値を用いて、漸化式にしたがって、収束解に近い複合材料の有効材料定数の値（ｎ＝２）を再帰的に求める。この値は、後述する収束解判定モジュール２６に送られ、判定の結果、解が収束しない場合、算出された材料定数の値を用いて、上記漸化式に従って、さらに収束解に近い複合材料の有効材料定数の値（ｎ＝３）を再帰的に求める。このようにして、方程式解法モジュール２４は、収束解判定モジュール２６において、収束解が見出されるまで、複合材料の有効材料定数の値を再帰的に求める処理を繰り返す。

また、第２の処理の場合、方程式解法モジュール２４は、下記に示す関数ｆ（ｘ）（ｘは求めようとする所定の材料の体積分率）を定め、この関数ｆ（ｘ）の上限値ｘ＝ｘ₁における値ｆ（ｘ₁）と下限値ｘ＝ｘ₂における値ｆ（ｘ₂）との間で積ｆ（ｘ₁）・ｆ（ｘ₂）を算出し、この結果を収束解判定モジュール２６に送る。

関数ｆ（ｘ）＝（実験により取得した複合材料の有効材料定数）−（再帰型非線形方程式（６）で求める体積分率ｘにおける複合材料の有効材料定数）

このような積ｆ（ｘ₁）・ｆ（ｘ₂）の算出は、収束解が見出されるまで、上限値ｘ＝ｘ₁と下限値ｘ＝ｘ₂における値を変更しながら繰り返し求められる。

収束解判定モジュール２６は、方程式解法モジュール２４の算出結果に基づいて収束解か否かを判定する部分である。
第１の処理の場合、繰返回数（ｎ＋１）で求めた複合材料の有効材料定数の値と、繰返回数ｎでもとめた複合材料の有効材料定数の値との差分の絶対値が予め設定された閾値を下回る場合、繰返回数（ｎ＋１）における複合材料の有効材料定数の値が収束解であるとされる。それ以外は、収束解は見出されないとして、繰返回数（ｎ＋１）で求めた値を用いて方程式解法モジュール２４で計算が行われるように、指示される。上述したように、再帰型非線型方程式は、係数等に具体的な数値が付与された方程式なので、この非線形方程式の微分した関数も既知である。したがって、収束解判定モジュール２６は、公知のニュートン・ラフソン法を用いて収束解を見出すことができる。
第２の処理の場合、ｆ（ｘ₁）・ｆ（ｘ₂）の値が、負か否かが調べられ、後述するように、２分法を用いて、上限値および下限値の値が再設定される。

結果処理モジュール２８は、収束解判定モジュール２６で得られた収束解を用いて、複合材料の有効ヤング率やの有効せん断剛性を算出し、あるいは、所定の材料の体積分率を取得する部分である。

出力装置３４は、結果処理モジュール２８で得られた各情報をプリント出力し、あるいは画面表示を行う部分である。
記憶装置３０には、予め設定された複合材料や複合材料に含まれる材料の材料定数等が記録保持されているデータベースを備えている。

なお、上述の方程式解法モジュール２４は、母相αに２つの異なる材料β，γを含む場合、母相αに材料β、γが分散した均一の材質の仮想複合材料を定めて、この仮想複合材料の有効材料定数を算出する。

このような処理は、複合材料の有効材料定数の算出を、以下に示す再帰型非線形方程式を用いて解析的に行うことにより実現される。以下、この再帰型非線形方程式について詳細に説明する。

図２（ａ）〜（ｃ）は、複合材料を説明する図である。
複合材料５０は、図２（ａ）に示すように、エポキシ樹脂（以降、エポキシという）Ｐからなる母相中に、ウレタンＵ及びフィラーＦを含む。
図２（ｂ）および（ｃ）は、複合材料５０に作用する応力およびひずみを示している。 以下で説明する応力、ひずみおよび材料定数は、それぞれ２階テンソル、２階テンソル、および４階テンソルであるが、わかり易く説明するためにスカラーの形式で説明する。
複合材料５０では、図２（ｂ）に示すように、ウレタンＵに作用する応力σ_Uと、フィラーＦに作用する応力σ_Fと、エポキシＰに作用する応力σ_Pとが加算されて、複合材料全体の平均応力σ_Tとなる。これに対応して、図２（ｃ）に示すように、ウレタンＵに作用するひずみε_Uと、フィラーＦに作用するひずみε_Fと、エポキシＰに作用するひずみε_Pとが加算されて、複合材料全体の平均ひずみε_Tとなる。平均応力σ_Tおよび平均ひずみε_Tは、下記式（１）、（２）で表される。平均応力及び平均ひずみは、いずれも体積分率を重み付け係数とした体積平均の応力及び体積平均のひずみを意味する。

一方、複合材料全体における平均応力σ_Tとひずみε_Tとの関係を規定する材料定数Ｃ^＊は、図３（ａ）に示すように、下記式（３）で表される。すなわち、材料定数Ｃ^＊を持つ均一な材質の仮想複合材料が想定される。

ここで、ひずみに注目し、ウレタンＵに作用するひずみε_Uと複合材料５０全体に作用するひずみε_Tとの間の関係を、下記式（４）で表す。同様に、フィラーＦに作用するひずみε_Fと複合材料５０全体に作用するひずみε_Tとの間の関係を、下記式（５）で表す。

式（４）は、ウレタンＵの平均ひずみであるε_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を用いて、式（１）〜（５）を整理すると、下記式（６）が導かれる。

式（６）は、左辺の材料定数Ｃ^＊を求める式であり、右辺にある比例定数Ａ^U，Ａ^Fはそれぞれ、Ｃ^U，Ｃ^＊およびＣ^F，Ｃ^＊に関して非線形な式となっている。このため、式（６）は、材料定数Ｃ^＊に関して、再帰型の非線形方程式となっている。

このような比例定数Ａ^U，Ａ^Fは、図４（ａ），（ｂ）に示すように、材料定数Ｃ^＊を持つ均一な材質の仮想複合材料中に球状粒子であるウレタンＵおよびフィラーＦが配置されるモデルを仮想し、この球状粒子を囲む周囲の材料定数を、算出すべき仮想複合材料の有効材料定数Ｃ^＊として定めることにより得られる。
すなわち、比例定数Ａ^U，Ａ^Fは、体積力が作用せず、無限遠で印加されるひずみε^∞ _ij（ｉ，ｊは、１〜３の自然数）が存在する条件下、公知のナビエの方程式を解くことにより得られる。なお、無限遠で印加されるひずみε^∞ _ijは下記式（７）のように、第１項の静水圧成分と第２項のせん断成分とに分解できる。このため、比例定数Ａ^U，Ａ^Fについても静水圧成分と第２項のせん断成分に分けて算出することができる。比例定数Ａ^U，Ａ^Fは、いずれも球状粒子の材料定数がＣ^UかＣ^Fの違いしかないので、以下の説明では、球状粒子の材料定数をＣ^U，Ｃ^Fの代表としてＣ^Bを用いて説明する。

以下、Ａ（Ｃ^B，Ｃ^＊）の第１項の静水圧成分と第２項のせん断成分の算出について説明する。

（静水圧成分に基づくＡ（Ｃ^B，Ｃ^＊）の算出）
式（７）中の第１項の静水圧成分に関するナビエの方程式を満足する変位ｕ_i（ｉは１〜３の自然数）は、下記式（８）のように表される。これは、右辺のテンソルの次元を変位ｕ_iの次元と対応させる必要性から一意的に定まる。
ここで、球状粒子における変位ｕ_iが原点で有限の値を持ち、ひずみが無限遠で式（７）中の第１項のひずみになる条件と、変位ｕ_iおよび球体粒子の表面における表面張力の連続性の条件とを与えることにより、球体粒子内部の変位ｕ_iと球体粒子外部の変位ｕ_iを求めることができる。具体的には、球体粒子内部の変位ｕ_iは下記式（９）で、球体粒子外部の変位ｕ_iは下記式（１０）で表される。

すなわち、母相に球状粒子がある空間内の各方向の変位ｕ_iは、球状粒子の内部において、球状粒子の中心を基準として定められる位置ｘ_i（i＝１、２または３）に比例する。変位ｕ_iは、さらに、球状粒子の外部において、球状粒子の中心を基準として定められる位置ｘ_i（i＝１、２または３）に比例する項と、球状粒子の中心からの距離の３乗に反比例し、位置ｘ_iに比例する項とを有するように定められる。これにより、比例定数Ａ（Ｃ^B，Ｃ^＊）の静水圧成分を求めることができる。

式（９）及び式（１０）から、球状粒子内部における応力σ_ijは下記式（１１）で表すことができ、球状粒子内部におけるひずみε_ijは、下記式（１２）で表すことができる。この式（１２）は、上述した式（５）の関係式の静水圧成分に対応するものである。したがって、Ａ（Ｃ^B，Ｃ^＊）の静水圧成分は、式（１３）のように表される。したがって、ウレタンＵの比例定数Ａ^Uは、式（１３）中のλ^Bおよびμ^Bをλ^Uおよびμ^Uに置き換えるとよい。同様に、フィラーＦの比例定数Ａ^Fは、式（１３）中のλ^Bおよびμ^Bをλ^Fおよびμ^Fに置き換えるとよい。

（せん断成分に基づくＡ（Ｃ^B，Ｃ^＊）の算出）
式（７）中の第２項のせん断成分に関するナビエの方程式を満足する変位ｕ_i（ｉは１〜３の自然数）は、球状粒子の内部及び外部を問わず、球状粒子の中心からの距離の０乗、距離の２乗、距離の−３乗および距離の５乗のそれぞれに比例する項の加算式で定められる。具体的には、変位ｕ_iを、｛ｆ（ｒ）・（ｘ_iｘ_j／ｒ²）＋ｇ（ｒ）・δij｝・ｘ_k（ｉ，ｊ，ｋ＝１〜３の整数）に比例するように定める。この式は、テンソルの階数をナビエの方程式上で対応させる必要性から一意的に定められたものである。ここで、ｘ_ｊ，ｘ_kは、３次元座標におけるｘ_iと異なる座標成分を表す。また、ｆ（ｒ）、ｇ（ｒ）は球状粒子の中心からの距離ｒのみによる多項式を表す。このとき、ナビエの方程式を満足するｆ（ｒ）、ｇ（ｒ）をｆ（ｒ）∝ｒ^ｌ、ｇ（ｒ）∝ｒ^ｌ（ｌは整数）としてｒの次数ｌを次元解析により求めることにより、ｒの次数ｌは０，２，−３，５となる。したがって、ナビエの方程式を満足する変位ｕ_iは、球状粒子の内部及び外部を問わず、球状粒子の中心からの距離の０乗、距離の２乗、距離の−３乗および距離の５乗のそれぞれに比例する項の加算式により表される。例えば、次数ｌ＝２のときの式は、下記式（１４）に示すように表される。勿論、式（１４）中のλ、μは、球状粒子の内部では、球状粒子の内部の材料定数が、球状粒子の外部では、球状粒子の外部の材料定数が用いられる。

すなわち、変位ｕ_iは、球状粒子の内部および外部において、球状粒子の中心からの距離の０乗、距離の２乗、距離の−３乗および距離の５乗のそれぞれに比例する項の加算式として表すことができる。これにより、再帰型非線形方程式を得るための比例定数Ａ（Ｃ^B，Ｃ^＊）のせん断成分を求めることができる。

具体的には、球状粒子の表面における変位ｕ_iおよび球体粒子の表面における表面張力の連続性の条件を与えて、下記式（１５）に示す球体粒子内部の応力が求められる。この式（１５）から、比例定数Ａ（Ｃ^B，Ｃ^＊）のせん断成分が下記式（１６）のように表される。したがって、ウレタンＵの比例定数Ａ^Uは、式（１６）中のμ^Bをμ^Uに置き換えるとよい。同様に、フィラーＦの比例定数Ａ^Fは、式（１６）中のμ^Bをμ^Fに置き換えるとよい。

以上より、ウレタンＵの比例定数Ａ^UおよびフィラーＦの比例定数Ａ^Fが定められ、これらの比例定数が、材料定数Ｃ^＊を求める式（６）に代入されて、再帰型非線形方程式が導かれる。この再帰型非線形方程式は、ウレタンＵおよびフィラーＦが球状粒子としてエポキシＰに分散された均質な材質の仮想複合材料を仮想し、この仮想複合材料中の球状粒子を囲む周囲の材料定数を、算出すべき仮想複合材料の有効材料定数として定めることにより実現したものである。
この仮想複合材料を利用した本発明の算出方法による複合材料の有効材料定数の計算結果は、後述するように、従来の長時間の計算を要する有限要素モデルを用いた場合に近く、精度が高いものとなっている。

このような再帰型非線形方程式を用いて複合材料の有効材料定数Ｃ^＊を算出する方法について説明する。

図５は、上記材料定数Ｃ^＊を算出するフローを説明する図である。
まず、条件設定モジュール２０にて、複合材料の有効材料定数Ｃ^＊を算出する第１の処理が設定される。この設定は、入力装置３２を介してオペレータから指示入力されて行われる。ここで、材料定数Ｃ^＊は、算出すべき未知数であるラメの定数λ_ef,μ_efを代表して表している。
次に、ウレタンＵ、フィラーＦおよびエポキシＰの材料定数Ｃ^U，Ｃ^F，Ｃ^P、すなわち、λ^P，μ^P，λ^U，μ^U，λ^F，μ^Fと、ウレタンＵおよびフィラーＦの体積分率ｖ_u，ｖ_Fの値が、記憶装置３０中のデータベースから呼び出されて取得される（ステップＳ１０）。材料定数と体積分率の値は、非線形方程式設定モジュール２２に送られて、上述の式（６）の各係数の値が定められ、材料定数Ｃ^＊のみの再帰型非線形方程式に修正される。すなわち、方程式が用意される。

次に、方程式解法モジュール２４にて、材料定数Ｃ^＊の未知数であるラメの定数λ_ef,μ_efの初期値が設定される（ステップＳ２０）。初期値を設定する理由は、式（６）が非線形であり、解が収束するまで繰り返すためである。初期値はどのような値であってもよいが、例えば、母相であるエポキシＰの材料定数Ｃ^PとフィラーＦの材料定数Ｃ^Fの、体積分率による加重平均値を選択するとよい。この初期値を、ラメの定数λ_ef,＝λ₁，μ_ef＝μ₁と定める。さらに、この初期値をそれぞれ、λ_n，μ_n（ｎ＝１）と表す。
次に、方程式解法モジュール２４にて、上述の式（６）を修正した下記式（１７）を定義する。具体的には、式（６）の左辺を右辺に移動し、その時の右辺をＦ（Ｃ^＊)として定義する。Ａ^ＵはＡ（Ｃ^Ｕ，Ｃ^＊）であり、Ａ^ＦはＡ（Ｃ^Ｆ，Ｃ^＊）である。

さらに、下記式（１８）にしたがって、Ｆ（Ｃ^＊)を用いてＣ⁽ⁿ⁺¹⁾が定義される。ここで、式（１８）は、Ｆ（Ｃ^＊)＝０となる収束解を求めるためにニュートン・ラフソン法を適用した式である。ここで、Ｃ⁽ⁿ⁾，Ｃ⁽ⁿ⁺¹⁾はラメの定数λ_n，μ_n，λ_n+1，μ_n+1を代表して表す符号である。式（１８）のＣ⁽ⁿ⁾に、ラメの定数λ_n，μ_nが与えられ、式（１８）に従ってラメの定数λ_n+1，μ_n+1が算出され取得される（ステップＳ３０）。ここで、式（１８）中のＦ´（Ｃ⁽ⁿ⁾）はヤコビのテンソルであり、Ｆ´（Ｃ⁽ⁿ⁾）^-1はＦ´（Ｃ⁽ⁿ⁾）の逆テンソルを表す。すなわち、Ｆ´（Ｃ(n)）^-1・Ｆ’（Ｃ⁽ⁿ⁾）＝Ｉ（Ｉは単位テンソル）である。

算出されたλ_n+1，μ_n+1と、この算出に用いられたλ_n，μ_nが収束解判定モジュール２６に送られて、それぞれの差分λ_n+1−λ_n，μ_n+1−μ_nの絶対値が求められる。この絶対値が予め設定された閾値ε₁，ε₂と比較される（ステップＳ４０）。比較の結果、差分λ_n+1−λ_n，μ_n+1−μ_nの絶対値がいずれも閾値ε₁，ε₂より小さい条件を満たす場合（Ｙｅｓの場合）、算出すべき仮想複合材料の有効材料定数Ｃ^＊であるラメの定数λ^*,μ^＊は、それぞれλ_n+1，μ_n+1に決定される（ステップＳ５０）。一方、比較の結果、上記条件を満たさない場合（Ｎｏの場合）、ｎ＋１をｎとし（ステップＳ６０）、さらに、式（１８）を用いてλ_n+1，μ_n+1を算出するためにステップＳ３０に戻る。こうして、ステップＳ４０における条件を満足するまで、ステップＳ３０，Ｓ４０，Ｓ６０が繰り返される。

次に、このように定められた複合材料の有効材料定数Ｃ^＊、具体的にはラメの定数λ^*,μ^＊を用いてヤング率Ｅ^*が算出される。ヤング率Ｅ^*は、下記式（１９）に従って算出され、プリンタやディスプレイ等の出力装置に出力される。

算出されたヤング率Ｅ^*は、後述するように、有限要素モデルを用いて数値計算をおこなって得られるヤング率Ｅとよく対応した値となっている。このため、本発明の複合材料の材料定数の算出方法は有効であることがわかる。また、本発明の複合材料の材料定数の算出方法は、式（６）に示される式を用いて解析的に算出できるので、有限要素モデルを用いた計算に比べて短時間に算出結果を得ることができ、効率的である。

本発明の複合材料の材料定数の算出方法で用いる式（６）は、解析式であるので、例えば、母相に含まれるフィラーＦ等の材料の体積分率を解くべき未知数とすることもできる。このとき、複合材料の有効材料定数は、予め実験等により求めておく必要がある。このような所定の材料の体積分率の算出は、上述した第２の処理に対応する。この第２の処理の具体的な方法を以下説明する。

図６は、第２の処理のフローを示す図である。エポキシＰに、フィラーＦとウレタンＵを含む複合材料において、フィラーＦの体積分率ｖ_Fを算出する場合を想定する。

まず、条件設定モジュール２０において、エポキシＰ、ウレタンＵ、フィラーＦの材料定数、具体的には、λ_P，μ_P，λ_U，μ_U，λ_F，μ_Fと、ウレタンＵの体積分率ｖ_Uが取得される（ステップＳ１００）。これらの材料定数等は、記憶装置３０中のデータベースから呼び出されて取得される。あるいは、入力装置３２から入力される。
次に、複合材料の力学的実験による計測結果より、複合材料全体の材料定数である有効材料、具体的にはλ^＊，μ^＊が特定され取得される（ステップＳ１１０）。計測結果は、予め記憶装置３０に記憶されており、条件設定モジュール２０に計測結果が呼び出されて、複合材料の有効材料定数が特定される。非線形方程式設定モジュール２２では、ステップＳ１００、Ｓ１１０にて取得された材料定数の各値が式（６）の各係数に代入されて、未知数が体積分率ｖ_Fの方程式が用意される。さらに、用意された式（６）を用いて、関数ｆ（ｘ）＝（実験により取得した複合材料の有効材料定数）−（式（６）で求める体積分率ｘにおける複合材料の有効材料定数）が定義される。

この後、フィラーＦの体積分率ｖ_Fの上限値ｘ₁及び下限値ｘ₂が設定される（ステップＳ１２０）。上限値および下限値は、後述する体積分率ｖ_Fの算出を、２分法を用いて行うために設定される。この設定は、入力装置３２を通して、オペレータの指示入力に基づいて設定されてもよいし、予め設定されたデフォルト値を用いてもよい。

次に、方程式解法モジュール２４において、上限値ｘ₁及び下限値ｘ₂が用いられて、ｆ（ｘ₁）・ｆ（ｘ₂）の値が負か否かが判定される（ステップＳ１３０）。ｆ（ｘ₁）・ｆ（ｘ₂）の値が正の場合、上限値ｘ₁及び下限値ｘ₂が変更される（ステップＳ１４０）。変更の方法は、特に制限されないが、ｘ₁及びｘ₂は体積分率であるので、上限値ｘ₁はより大きく、下限値ｘ₂はより小さくなるように設定されるとよい。

ステップＳ２３０における判定で肯定される場合（Ｙｅｓの場合）、公知の２分法に基づいて、以下の処理が行われる。
すなわち、まず、ｘ₃＝（ｘ₁＋ｘ₂）／２が求められ（ステップＳ１５０）、さらに、ｆ（ｘ₁）・ｆ（ｘ₃）の値が負か否かが判定される（ステップＳ１６０）。判定の結果、ｆ（ｘ₁）・ｆ（ｘ₃）が負である場合（Ｙｅｓの場合）、下限値ｘ₂がｘ₃に変更され（ステップＳ１７０）、ｆ（ｘ₁）・ｆ（ｘ₃）が負でない場合（Ｎｏの場合）、上限値ｘ₁がｘ₃に変更される（ステップＳ１８０）。
この後、上限値ｘ₁と下限値ｘ₂との差分の絶対値が予め設定された閾値ε₅より小さいか否かが判定される（ステップＳ１９０）。絶対値が閾値ε₅より小さい場合（Ｙｅｓの場合）、上限値ｘ₁または下限値ｘ₂が算出すべき体積分率ｖ_Fとして定められる（ステップＳ２００）。ステップＳ１９０において、絶対値が閾値ε₅より小さくない場合（Ｎｏの場合）、ステップＳ１５０に戻る。こうして、ステップＳ１９０の判定で、絶対値が閾値ε₅より小さくなるまで、ステップＳ１５０、Ｓ１６０、Ｓ１７０、Ｓ１８０を繰り返す。

以上の方法により、フィラーＦの体積分率ｖ_Fは算出される。算出したフィラーＦの体積分率ｖ_Fは、出力装置３４に出力される。
なお、算出すべき体積分率は、フィラーＦを対象とする他、ウレタンＵを対象とすることもできる。また、算出すべき体積分率は、複数の種類の材料を対象としてもよい。
複合材料の材料定数の算出方法および体積分率の算出方法では、母相に２つの材料が分散した複合材料を用いて説明したが、３つ以上であってもよい。母相は、エポキシに限られず、金属やセラミック等の無機材料であってもよい。また、母相に含まれる材料は、母相を強化する強化材料に限られず、上述したようにウレタン等の軟質材料であってもよい。

また、複合材料の材料定数の算出方法および体積分率の算出方法は、コンピュータ上でプログラムを実行することにより実現することもできる。このプログラムはコンピュータで読み取り可能な記録メディアに記録されたものである。この記録メディアは、通信回線を利用してダウンロードされたプログラムを記録したものも含まれる。

次に、上述の複合材料の有効材料定数の算出方法の有効性を説明する。
複合材料は、以下に示すものが用いられた。
（１）エポキシ、ウレタン：６０ｇ
ウレタン変性エポキシ： ３０ｇ
ＮＢＲ変性エポキシ ： ２０ｇ
ビスフェノールＡ型液状エポキシ：１０ｇ
（２）フィラー：３５ｇ
シリカ ：１５ｇ
炭酸カルシウム ：１０ｇ
酸化カルシウム ：１０ｇ
（３）硬化系：５ｇ
ＤＩＣＹ（ジシアンジアミド） ：５ｇ

複合材料におけるヤング率及びポアソン比と、体積分率（％）は以下の通りである。なお、ヤング率およびポアソン比から、公知の式により、ラメの定数に変換して用いた。

本発明の複合材料の材料定数の算出方法を用いた複合材料の有効材料定数として得られるヤング率は２．１２６（ＧＰａ）であった。一方、有限要素モデルを用いた複合材料の解析により得られたヤング率は２．００７（ＧＰａ）であった。
これより、本発明の方法により得られるヤング率は、有限要素モデルを用いて得られるヤング率および実験により得られたヤング率と極めて近い値を示し、本発明の方法が有効であることがわかる。

なお、有限要素モデルを用いた材料定数の算出は、図７に示すような、エポキシＰを母相とし、このエポキシＰ中に球状のウレタンＵが体積分率３２％含まれ、フィラーＦが体積分率１４％含まれて構成される複合材料を想定した有限要素モデルを作成して、複合材料の有効材料定数を算出した。図中上端部に引っ張り変位を与えて応力ーひずみの初期傾きを調べてヤング率とした。なお、作成した有限要素モデルは、要素数６５５３６個、節点数６６０４９個の２次元平面応力要素モデルである。解析は、非線形有限要素法ソフトウェアＡＢＡＱＵＳ（商品名）を用いた。

以上、本発明の複合材料の材料定数の算出方法、複合材料中の材料の体積分率の算出方法および記録メディアについて詳細に説明したが、本発明は上記実施形態に限定されず、本発明の主旨を逸脱しない範囲において、種々の改良や変更をしてもよいのはもちろんである。

１０ 計算装置
１２ ＣＰＵ
１４ ＲＯＭ
１６ ＲＡＭ
１８ 入出力ポート
２０ 条件設定モジュール
２２ 非線形方程式設定モジュール
２４ 方程式解法モジュール
２６ 収束解判定モジュール
２８ 結果処理モジュール
３０ 記憶装置
３２ 入力装置

Claims (12)

1. 力学的材料定数が既知の複数の種類の材料を、力学的材料定数が既知の母相中に含む複合材料に関して、複合材料の力学的有効材料定数をコンピュータが算出する方法であって、
   前記母相中に、前記複数の種類の材料のそれぞれが既知の体積分率で球状粒子として分散した仮想複合材料を、前記複合材料として定めることにより、この仮想複合材料の有効材料定数を未知数とする非線形方程式を用意するステップと、
   前記用意した非線形方程式を解いて、前記仮想複合材料の有効材料定数を、前記複合材料の有効材料定数として求めるステップと、を有し、
   前記非線形方程式は、前記仮想複合材料中の前記複数の種類の材料のそれぞれに対応する前記球状粒子を囲む周囲の材料定数を、算出すべき前記仮想複合材料の有効材料定数として定めることにより得られる再帰型非線型方程式である複合材料の有効材料定数の算出方法。
2. 前記母相の材料定数をＣ^Aとし、前記複数の種類の材料の材料定数をＣ^B1，Ｃ^B2，………Ｃ^Bn(ｎは２以上の整数)とし、前記複数の材料の体積分率をそれぞれｖ_B1，ｖ_B2，………ｖ_Bnとし、前記仮想複合材料の有効材料定数をＣ^*としたとき、
   前記非線形方程式は下記式で表され、
   下記式中のＡ（Ｃ^B1，Ｃ^＊），Ａ（Ｃ^B2，Ｃ^＊），………，Ａ（Ｃ^Bn，Ｃ^＊）は、前記仮想複合材料中の前記球状粒子を囲む周囲の材料定数を、算出すべき前記仮想複合材料の有効材料定数として定めることにより得られ、前記材料定数Ｃ^＊に関して非線形な式で定まる比例定数である、請求項１に記載の算出方法。
   Ｃ^* ＝ Ｃ^A＋ｖ_B1・（Ｃ^B1−Ｃ^A）・Ａ（Ｃ^B1，Ｃ^＊）
   ＋ｖ_B2・（Ｃ^B2−Ｃ^A）・Ａ（Ｃ^B2，Ｃ^＊）
   ………＋ｖ_Bn・（Ｃ^Bn−Ｃ^A）・Ａ（Ｃ^Bn，Ｃ^＊）
3. 前記仮想複合材料の応力およびひずみを、せん断成分と静水圧成分に分解し、前記仮想複合材料中の、前記静水圧成分に対応する変位をｕ_i（i＝１、２または３）としたとき、
   前記変位ｕ_iは、前記球状粒子の内部において、前記球状粒子の中心を基準とする位置ｘ_i（i＝１、２または３）に比例し、
   前記変位ｕ_iは、前記球状粒子の外部において、前記球状粒子の中心を基準とする位置ｘ_i（i＝１、２または３）に比例する項と、前記球状粒子の中心からの距離の３乗に反比例し、前記位置ｘ_i（i＝１、２または３）に比例する項を有するように定めることにより、前記非線形方程式が用意される、請求項１または２に記載の算出方法。
4. 前記仮想複合材料の応力およびひずみを、せん断成分と静水圧成分に分解し、前記仮想複合材料中の、前記せん断成分に対応する変位をｕ_i（i＝１、２または３）としたとき、
   前記変位ｕ_iを、前記球状粒子の内部および外部において、前記球状粒子の中心からの距離の０乗、距離の２乗、距離の−３乗および距離の５乗のそれぞれに比例する項の加算式で定めることにより、前記非線形方程式が用意される、請求項１または２に記載の算出方法。
5. 前記非線形方程式は、ニュートン・ラフソン法を用いて、解が収束するように解かれる請求項１〜４のいずれか１項に記載の算出方法。
6. 力学的材料定数が既知の母相中に、力学的材料定数が既知の複数の種類の材料を含む複合材料に関して、いずれか１種類の材料の体積分率をコンピュータが算出する方法であって、
   前記複合材料の有効材料定数を実験結果から特定するステップと、
   前記母相中に、前記複数の種類の材料のそれぞれが球状粒子として分散した仮想複合材料を、前記複合材料として定めることにより、この仮想複合材料中のいずれか１種類の材料の体積分率を未知数とする非線形方程式を用意するステップと、
   前記用意した非線形方程式を解いて、前記材料の体積分率を求めるステップと、を有し、
   前記非線形方程式は、前記仮想複合材料中の前記複数の種類の材料のそれぞれの前記球状粒子を囲む周囲の材料定数を、前記特定された前記複合材料の有効材料定数として定めることにより得られる再帰型非線型方程式である、算出方法。
7. 前記母相の材料定数をＣ^Aとし、前記複数の種類の材料の材料定数をＣ^B1，Ｃ^B2，………Ｃ^Bn(ｎは２以上の整数)とし、前記仮想複合材料の有効材料定数をＣ^*とし、前記複数の材料の体積分率をそれぞれｖ_B1，ｖ_B2，………ｖ_Bnとしたとき、
   前記非線形方程式は下記式で表され、かつ、Ａ（Ｃ^B1，Ｃ^＊），Ａ（Ｃ^B2，Ｃ^＊），………，Ａ（Ｃ^Bn，Ｃ^＊）は、前記仮想複合材料中の前記球状粒子を囲む周囲の材料定数を、算出すべき前記仮想複合材料の有効材料定数として定めることにより得られ、前記材料定数Ｃ^＊に関して非線形な式で定まる比例定数である、請求項６に記載の算出方法。
   Ｃ^* ＝ Ｃ^A＋ｖ_B1・（Ｃ^B1−Ｃ^A）・Ａ（Ｃ^B1，Ｃ^＊）
   ＋ｖ_B2・（Ｃ^B2−Ｃ^A）・Ａ（Ｃ^B2，Ｃ^＊）
   ………＋ｖ_Bn・（Ｃ^Bn−Ｃ^A）・Ａ（Ｃ^Bn，Ｃ^＊）
8. 前記仮想複合材料の応力およびひずみを、せん断成分と静水圧成分に分解し、前記仮想複合材料中の、前記静水圧成分に対応する変位をｕ_i（i＝１、２または３）としたとき、
   前記変位ｕ_iは、前記球状粒子の内部において、前記球状粒子の中心を基準とする位置ｘ_i（i＝１、２または３）に比例し、
   前記変位ｕ_iは、前記球状粒子の外部において、前記球状粒子の中心を基準とする位置ｘ_i（i＝１、２または３）に比例する項と、前記球状粒子の中心からの距離の３乗に反比例し、前記位置ｘ_i（i＝１、２または３）に比例する項を有するように定めることにより、前記非線形方程式が用意される、請求項６または７に記載の算出方法。
9. 前記仮想複合材料の応力およびひずみを、せん断成分と静水圧成分に分解し、前記仮想複合材料中の、前記せん断成分に対応する変位をｕ_i（i＝１、２または３）としたとき、
   前記変位ｕ_iを、前記球状粒子の内部および外部において、前記球状粒子の中心からの距離の０乗、距離の２乗、距離の−３乗および距離の５乗のそれぞれに比例する項の加算式で定めることにより、前記非線形方程式が用意される、請求項６または７に記載の算出方法。
10. 前記非線形方程式は、２分法を用いて解かれる請求項６〜９のいずれか１項に記載の算出方法。
11. 力学的材料定数が既知の母相中に、力学的材料定数が既知の複数の種類の材料を含む複合材料に関して、複合材料の力学的有効材料定数を算出する、コンピュータの実行可能なプログラムを記録した記録メディアであって、
    前記プログラムは、
    前記母相中に、前記複数の種類の材料のそれぞれが既知の体積分率で球状粒子として分散した仮想複合材料を、前記複合材料として定めることにより、この仮想複合材料の有効材料定数を未知数とする非線形方程式を用意する手順と、
    前記用意した非線形方程式を解いて、前記仮想複合材料の有効材料定数を、前記複合材料の有効材料定数として求める手順とを、コンピュータに実行させ、
    前記非線形方程式は、前記仮想複合材料中の前記複数の種類の材料のそれぞれに対応する前記球状粒子を囲む周囲の材料定数を、算出すべき前記仮想複合材料の有効材料定数として定めることにより得られる再帰型非線型方程式である、記録メディア。
12. 力学的材料定数が既知の母相中に、力学的材料定数が既知の複数の種類の材料を含む複合材料に関して、いずれか１種類の材料の体積分率を算出する、コンピュータの実行可能なプログラムを記録した記録メディアであって、
    前記プログラムは、
    前記複合材料の有効材料定数を実験結果から特定する手順と、
    前記母相中に、前記複数の種類の材料のそれぞれが球状粒子として分散した仮想複合材料を、前記複合材料として定めることにより、この仮想複合材料中のいずれか１種類の材料の体積分率を未知数とする非線形方程式を用意する手順と、
    前記用意した非線形方程式を解いて、前記材料の体積分率を求める手順とを、コンピュータに実行させ、
    前記非線形方程式は、前記仮想複合材料中の前記複数の種類の材料のそれぞれの前記球状粒子を囲む周囲の材料定数を、前記特定された前記複合材料の有効材料定数として定めることにより得られる再帰型非線型方程式である、記録メディア。