(57)【要約】      （修正有）【課題】少ない要素数で圧電磁器の解析を行うことができるとともに、未分極状態から分極状態への遷移過程を表現することで非線形圧電現象を考慮できる圧電磁器の物理量算出方法を提供する。【解決手段】圧電磁器の物性値および自発ひずみを用いて、一つの有限要素を構成する複数の結晶方位のセルのドメインウォールが変化する駆動力の式および非線形項を取り入れた圧電基本式とから求めた仮想のバタフライカーブ、仮想のＤ−Ｅヒステリシスカーブが、実測のカーブに近似するように、駆動力の式のパラメータを求め、多数の有限要素について計算実行し、圧電磁器の物理量を求める方法であって、結晶方位が複数のセルから構成され、該複数のセルについてそれぞれ駆動力の式の値を求め、該駆動力の式の値がある一定値を超えたセルについてはドメインウォールを変化させて、仮想のバタフライカーブ、仮想のＤ−Ｅヒステリシスカーブを求める。【選択図】図１
2010258188.tif 000002
【特許請求の範囲】
【請求項１】
圧電磁器が有する弾性スティフネスマトリクス、圧電ｅ定数マトリクス、誘電率マトリクスからなる物性値を得るとともに、自発ひずみ、バタフライカーブおよびＤ−Ｅヒステリシスカーブを実測して得る第１ステップと、
前記圧電磁器の物性値および自発ひずみを用いて、一つの有限要素を構成する複数の結晶方位のセルのドメインウォールが変化する駆動力の式と非線形項（自発分極、自発ひずみ）を取り入れた圧電基本式とから求めた仮想のバタフライカーブおよび仮想のＤ−Ｅヒステリシスカーブが、前記圧電磁器から得られた実測の前記バタフライカーブおよび前記Ｄ−Ｅヒステリシスカーブに近似するように、前記駆動力の式のパラメータを求める第２ステップと、
前記パラメータを用いた前記駆動力の式と前記非線形項（自発分極、自発ひずみ）を取り入れた圧電基本式とを用いて、多数の有限要素について、該多数の有限要素の少なくとも一つに電場および荷重の少なくとも一種を印加し、もしくは前記多数の有限要素の少なくとも一つに変位を与えて有限要素法で計算を実行し、前記圧電磁器の応力、分極およびひずみのうち少なくとも一種を求める第３ステップとを具備する圧電磁器の物理量算出方法であって、
前記第２ステップ中において、前記結晶方位が複数のセルから構成され、該複数のセルについてそれぞれ前記駆動力の式の値を求め、該駆動力の式の値がある一定値を超えたセルについてはドメインウォールを変化させて、前記仮想のバタフライカーブおよび前記仮想のＤ−Ｅヒステリシスカーブを求めることを特徴とする圧電磁器の物理量算出方法。
【請求項２】
前記第２ステップ中において、前記駆動力の式の値がある一定値を超えた場合には、前記セルのドメインウォールの方向変化または前記セルのドメインウォールの平行移動が発生したとして、前記仮想のバタフライカーブおよび前記仮想のＤ−Ｅヒステリシスカーブを求めることを特徴とする請求項１記載の圧電磁器の物理量算出方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、圧電磁器の物理量算出方法に関し、特に、圧電磁器に大きな電場や応力を印加した際に生じる非線形圧電現象を考慮した圧電磁器の物理量算出方法に関するものである。
【背景技術】
【０００２】
圧電磁器は強誘電体ドメインが動くことにより、電場や応力に対して非線形な応答を示すこと（非線形圧電現象）が知られている。非線形圧電現象は圧電磁器に大きな電場を印加した場合に顕著であり、電場に対する電束密度、ひずみの非線形な変化を描いたカーブはそれぞれＤ−Ｅヒステリシスカーブ、バタフライカーブと呼ばれる。この非線形圧電現象が生じる実際のデバイスとしては、大きな変位を生じさせるために大きな電場や応力を印加する圧電アクチュエータが知られている。
【０００３】
非線形圧電現象は圧電磁器部品の機能を高めるためには有用であるが、一方で圧電磁器部品の機能劣化やクラックの原因となっている。従って大きな電場や応力を印加する必要のある圧電デバイスでは、非線形圧電現象を考慮した高度な設計が必要である。
【０００４】
一般に大きな変位を必要とする圧電デバイスでは、ＰｂＺｒＯ３とＰｂＴｉＯ３の固溶体であるＰＺＴが用いられる。このＰＺＴの圧電定数は、結晶系が菱面晶相（ＰｂＺｒＯ３リッチ側）と正方晶相（ＰｂＴｉＯ３リッチ側）との境界であるモルフォトロピック相境界（ＭＰＢ）においてピークを示すことが知られている。ＭＰＢでは大きな圧電定数が得られる一方、強誘電体ドメインが動きやすく、非線形圧電現象も顕著である。
【０００５】
大きな電場や応力を印加する圧電デバイスでは非線形圧電現象が大きな影響を及ぼすため、この現象を考慮することのできる設計技術、シミュレーション技術が必要である。しかし非線形圧電現象を考慮できるプログラムは現在市販されておらず、研究者によっていくつかの計算手法が提案されている状況である。例えば、１つの要素が単独の強誘電体ドメインであり、その分極方向を電場や応力に応じて９０°もしくは１８０°変化させる有限要素法の計算手法が提案されている（非特許文献１参照）。また、強誘電体ドメインの動きを考慮せず、現象論的なパラメータを導入することで非線形圧電現象を表現する有限要素法の計算手法が提案されている（非特許文献２参照）。
【先行技術文献】
【非特許文献】
【０００６】
【非特許文献１】Mechanics of Materials, 36, pp.959-973, 2004
【非特許文献２】International Journal of Solid and Structure, 43, pp.7197-7222, 2006
【発明の概要】
【発明が解決しようとする課題】
【０００７】
しかしながら、非特許文献１に開示された計算手法では、１つの有限要素をシングルドメインとして扱っているため、Ｄ−Ｅヒステリシスカーブやバタフライカーブを表現するのに１，０００個程度の有限要素を必要としている。従ってこの手法では個々の要素の表現力が不足しており、実際の圧電デバイスを計算するには非現実的な要素数を必要とするという問題点があった。
【０００８】
また非特許文献２に開示された計算手法では、強誘電体ドメインの動きを考慮せず、現象論的なパラメータを導入することで１つの有限要素でＤ−Ｅヒステリシスカーブ、バタフライカーブを表現している。しかしこの手法では、電場を印加して未分極状態から分極状態となる遷移過程は表現できておらず、Ｄ−Ｅヒステリシスカーブ、バタフライカーブでは、未分極状態から分極状態への遷移過程では電束密度、ひずみが０となっている。
【０００９】
すなわち、実際の圧電アクチュエータでは電極端付近に疲労やクラック等の問題が生じることが多く、この領域では種々の程度に分極された領域が存在しており、この領域の応力を正しく求めるには周辺部の未分極状態の要素から電極端付近の大きく分極された要素までを計算する必要があるため、この手法では電極端先端の応力状態を上手く計算できないという問題点があった。
【００１０】
本発明は、少ない要素数で圧電磁器の解析を行うことができるとともに、未分極状態から分極状態への遷移過程を表現することで非線形圧電現象を考慮できる圧電磁器の物理量算出方法を提供することを目的とする。
【課題を解決するための手段】
【００１１】
本発明の圧電磁器の物理量算出方法は、圧電磁器が有する弾性スティフネスマトリクス、圧電ｅ定数マトリクス、誘電率マトリクスからなる物性値を得るとともに、自発ひずみ、バタフライカーブおよびＤ−Ｅヒステリシスカーブを実測して得る第１ステップと、
前記圧電磁器の物性値および自発ひずみを用いて、一つの有限要素を構成する複数の結晶方位のセルのドメインウォールが変化する駆動力の式と非線形項（自発分極、自発ひずみ）を取り入れた圧電基本式とから求めた仮想のバタフライカーブおよび仮想のＤ−Ｅヒステリシスカーブが、前記圧電磁器から得られた実測の前記バタフライカーブおよび前記Ｄ−Ｅヒステリシスカーブに近似するように、前記駆動力の式のパラメータを求める第２ステップと、
前記パラメータを用いた前記駆動力の式と前記非線形項（自発分極、自発ひずみ）を取り入れた圧電基本式とを用いて、多数の有限要素について、該多数の有限要素の少なくとも一つに電場および荷重の少なくとも一種を印加し、もしくは前記多数の有限要素の少なくとも一つに変位を与えて有限要素法で計算を実行し、前記圧電磁器の応力、分極およびひずみのうち少なくとも一種を求める第３ステップとを具備する圧電磁器の物理量算出方法であって、
前記第２ステップ中において、前記結晶方位が複数のセルから構成され、該複数のセルについてそれぞれ前記駆動力の式の値を求め、該駆動力の式の値がある一定値を超えたセルについてはドメインウォールを変化させて、前記仮想のバタフライカーブおよび前記仮想のＤ−Ｅヒステリシスカーブを求めることを特徴とする。
【００１２】
本発明の圧電磁器の物理量算出方法では、１つの有限要素が仮想的に複数の結晶方位の領域に分けられ、1つの有限要素の物理量は複数の結晶方位の寄与を合計したものであり、複数の結晶方位はそれぞれ異なる自発分極、自発ひずみを持つ複数の仮想的なセルより構成され、各セルの体積比率は外部電場、外部応力及びドメインウォールの駆動力によって非線形に変化し、結果として1つの有限要素で非線形圧電現象を表現することができる。
【００１３】
この計算手法では、１つの有限要素で多数のドメインの動きを計算しており、ある大きさの電場を印加した場合に一部のドメインは動くがその他は動かない、といったことが可能である。このため、１つの有限要素でＤ−Ｅヒステリシスカーブ、バタフライカーブや未分極から分極状態への遷移過程を表現することができる。
【００１４】
また、本発明の圧電磁器の物理量算出方法は、前記第２ステップ中において、前記駆動力の式の値がある一定値を超えた場合には、前記セルのドメインウォールの方向変化または前記セルのドメインウォールの平行移動が発生したとして、前記仮想のバタフライカーブおよび前記仮想のＤ−Ｅヒステリシスカーブを求めることを特徴とする。これにより、非線形圧電現象を精度良く表現することができる。
【発明の効果】
【００１５】
本発明によれば、少ない要素数で圧電磁器の物理量を、非線形圧電現象を考慮して精度良く算出できるとともに、未分極状態から分極状態への遷移過程を表現することができる。
【図面の簡単な説明】
【００１６】
【図１】本発明の圧電磁器の物理量算出方法のフローチャートである。
【図２】本発明の圧電磁器の物理量算出方法の第２ステップを詳細に説明するためのフローチャートである。
【図３】強誘電体ドメイン構造の説明図であり、（ａ）はＰＺＴのＴＥＭ（透過電子顕微鏡）写真、（ｂ）はドメイン構造の模式図、（ｃ）はＰＺＴ単位格子の模式図である。
【図４】本発明に係るドメインウォールの２つの動き、方向変化と平行移動の模式図である。
【図５】強誘電体ドメイン構造の説明図であり、ａ軸とｃ軸を区別しない場合のドメインウォールのとり得る方向の模式図である。
【図６】ＰＺＴのＤ−Ｅヒステリシスカーブ、バタフライカーブについて実測と作成した有限要素法プログラムでの計算結果とを比較した図である。
【発明を実施するための形態】
【００１７】
本発明の圧電磁器の物理量算出方法について、図１を用いて説明する。
【００１８】
第１ステップ
先ず、圧電磁器の各種物性値および駆動力のパラメータを求める。
【００１９】
物理量を算出しようとする圧電磁器の各種物性値を求める。求める物性値は、弾性スティフネスマトリクス〔ｃ〕、圧電ｅ定数マトリクス〔ｅ〕、誘電率マトリクス〔ε〕であり、ＥＭＡＳ試験（ＥＭＡＳ−６１００、圧電セラミック振動子の電気的試験方法、日本電子材料工業会標準規格）により求める。
【００２０】
また、圧電磁器のドメインの自発ひずみの大きさγＳＰを、Ｘ線回折測定にて求める。
【００２１】
さらに、圧電磁器のＤ−Ｅヒステリシスカーブ、バタフライカーブを測定する。
【００２２】
第２ステップ
次に、1つの有限要素（要素ということもある）で仮想的なＤ−Ｅヒステリシスカーブ、バタフライカーブを計算し、その仮想のＤ−Ｅヒステリシスカーブ、バタフライカーブが、実測のＤ−Ｅヒステリシスカーブ、バタフライカーブに近似するようにパラメータをフィッティングする。 
【００２３】
具体的には、図２に示すように、初期では要素は未分極とし、要素に与える電場を変化させながら要素に生じる電束密度、ひずみを計算する。
【００２４】
図2のフローチャートについて順を追って説明する。まず未分極の要素にある小さな電場を与え、下記の非線形項（自発分極、自発ひずみ）を取り入れた圧電基本式（ｅ形式）より要素の電場、応力、ひずみ、電束密度を計算する。
【００２５】
【数１】
2010258188.tif 000003
【００２６】
未分極の状態では非線形項である自発分極、自発ひずみ、および圧電ｅ定数マトリクスは共に０である。次に得られた電場、応力よりドメインウォールの方向変化、平行移動の駆動力Ｆｄ、Ｆｐを計算し、ドメインウォールの状態を求める（ドメインウォールの計算については後述の「第２ステップの詳細説明」を参照）。ドメインウォールの状態より要素の自発分極、自発ひずみ、および圧電ｅ定数マトリクスを計算することができる。
【００２７】
電場が小さい場合にはドメインウォールの状態は変化しないため反復計算は不要であり、次の電場の計算に移る。電場が大きくなってくるとドメインウォールの状態が大きく変化するため、収束計算が必要になり、1つの電場に対する収束計算の回数が多くなる。電場を変化させながら上述の計算を繰り返し、各電場での要素の｛Ｅ｝、｛Ｔ｝、｛D｝、｛Ｓ｝、自発分極、自発ひずみ等の物理量を得る。
【００２８】
なお、Ｄ−Ｅヒステリシスカーブでの縦軸は要素としての電束密度であり、線形項｛D｝と自発分極の和である。またバタフライカーブの縦軸は同じく要素としてのひずみであり、線形項｛Ｓ｝と自発ひずみの和である。この結果、仮想的なＤ−Ｅヒステリシスカーブ、バタフライカーブが得られる。次に得られた仮想的なＤ−Ｅヒステリシスカーブ、バタフライカーブと実測したＤ−Ｅヒステリシスカーブ、バタフライカーブとを比較し、実測したＤ−Ｅヒステリシスカーブ、バタフライカーブに近い結果が得られるようパラメータを変更しながら仮想的なＤ−Ｅヒステリシスカーブ、バタフライカーブを計算する。この繰り返しにより、実測したカーブに近い仮想的なＤ−Ｅヒステリシスカーブ、バタフライカーブを計算するためのパラメータを算出する。
【００２９】
第３ステップ
第１ステップ、第２ステップで得られた各種物性値及びパラメータを用いて多数の有限要素にて計算を実行し、応力、分極およびひずみの分布を求める。
【００３０】
第２ステップについて詳細に説明する。
【００３１】
図３は本発明に係る強誘電体ドメイン構造の説明図であり、（ａ）はＰＺＴの透過電子顕微鏡（ＴＥＭ）写真、（ｂ）はドメイン構造の模式図、（ｃ）はＰＺＴ単位格子の模式図である。ここでの説明はＰＺＴが正方晶であるとして行う。
【００３２】
正方晶のＰＺＴ単位格子には（Ｚｒ，Ｔｉ）原子の位置によりａ軸とｃ軸があり、自発分極の方向がｃ軸方向である。また自発ひずみによりｃ軸方向はａ軸方向と比べて１％程度長くなっている。
【００３３】
また、図３（ｂ）に示すように、単位格子の分極方向が揃っている領域をドメインと呼び、ドメインはドメインウォールにより区切られている。このドメインウォールを挟んでドメインの分極方向が９０°異なる場合には、このドメインウォールのことを９０°ドメインウォールという。正方晶ＰＺＴには９０°ドメインウォールと１８０°ドメインウォールが存在することが知られているが、ここでは１８０°ドメインについては考慮しない。
【００３４】
ＰＺＴ多結晶体は図３（ａ）に示すように個々の結晶粒がドメインウォールの方向が揃った複数の領域に区切られており、このドメインウォールの方向が揃った領域のことを、本願ではセルという。一般的にＰＺＴは多結晶体で使用され、結晶粒子は圧電デバイスのサイズと比べて非常に小さいため、圧電デバイスの解析を行う際には多数の結晶粒子を含む領域を計算する必要がある。
【００３５】
本発明の圧電磁器の物理量算出方法では、ドメインの動きとして、図４に示すように９０°ドメインウォールの法線ベクトルの変化（以下ドメインウォールの方向変化ということがある）とドメインウォールの平行移動を考える。なお、ここでのドメインウォールの法線ベクトルとは、ドメインウォールの２つの法線ベクトルのうち、ドメインウォールが形成されている領域の平均的な分極方向に近い方を指す。電場を印加した場合には、（図４（ａ））の状態から、セルの平均的な分極方向が印加電場の方向に近くなるようにドメインウォールの方向が変化することと（図４（ｂ））、印加電場に近い自発分極を持つドメインの領域が広がるようにドメインウォールが平行移動すること（図４（ｃ））の２通りの動きが生じる。
【００３６】
なお、図５に示すように、正方晶のａ軸とｃ軸を区別しない場合には、１つの結晶方位についてドメインの自発分極の方向は６通りである。９０°ドメインウォールの法線ベクトルは、ドメインウォールが区切る２つのドメインの体積が等しい場合のセルの平均的な分極方向に等しく、そのとり得る方向は１２通りである。
【００３７】
例えば、ａ軸、ｃ軸がｘ、ｙ、ｚ軸に平行である方位を持つ結晶のドメインウォールの法線ベクトルは、具体的には、以下で表される。
【００３８】
【数２】
2010258188.tif 000004
【００３９】
ドメインウォールの方向はどの方向にでも自由に変化できるのではなく、とり得る方向の中で電場や応力に応じて変化する。正方晶ではドメインウォールの変化する角度は６０°、９０°、１２０°、１８０°の４通りである。
【００４０】
また本発明の圧電磁器の物理量算出方法では、１つの要素で複数の結晶方位を考慮し、個々の結晶方位での計算結果の寄与を合計したものが要素の計算結果となる。またそれぞれの結晶方位は上述のように６通りの自発分極を持つドメイン、及び９０°ドメインウォールを形成する２つのドメインの組み合わせにより１２通りの平均的な分極方向を持つセルより構成されている。未分極状態では個々の結晶方位について１２のセルの体積比率は等しく、また個々のセルについてドメインウォールが区切る２つのドメインの体積比率は等しい。未分極状態に電場を印加すると、一部のセルのドメインウォールの方向変化によりセルの体積比率が変化し、またドメインウォールの平行移動によりドメインウォールが区切る２つのドメインの体積比率が変化する。
【００４１】
さらに本発明によれば、あるセルのドメインウォールの法線ベクトルが、例えば、図４（ａ）に示すｉ方向から図４（ｂ）に示すｊ方向へと変化する駆動力Ｆｄは以下の式で計算される。
【００４２】
【数３】
2010258188.tif 000005
【００４３】
Ｆｄの値がある閾値を超えた場合には、ドメインウォールの方向変化が発生し、ドメインウォールの法線ベクトルがｉ方向のセルの体積比率は０になり、その分ｊ方向の法線ベクトルを持つセルの体積比率が増加する。ここで、ドメインウォールの方向変化の閾値は、４通りのドメインウォールの方向変化に対応し４つ存在する。また、Ｆｄの計算では
【００４４】
【数４】
2010258188.tif 000006
【００４５】
の値は常にドメインウォールが２つのドメインを等しい体積に区切って計算することが望ましい。
【００４６】
また本発明によれば、ｆ方向とｇ方向の自発分極を持つドメインを区分する９０°ドメインウォールが、ｇ方向の自発分極を持つドメインの領域が広がるように平行移動する場合の駆動力Ｆｐは以下の式で計算される。
【００４７】
【数５】
2010258188.tif 000007
【００４８】
Ｆｐがある閾値を超えた場合には、ドメインウォールの平行移動が発生し、ｆ方向の自発分極を持つドメインの体積比率が減少し、その分ｇ方向の自発分極を持つドメインの体積比率が増加する。
【００４９】
これらのドメインウォールの動きは図２のフローチャートにて前述の（１）、（２）式の計算後に計算され、各セルの自発分極、自発ひずみおよび圧電ｅ定数マトリクスの寄与を全て計算して体積比率を乗じて和を取ることで要素の自発分極、自発ひずみおよび圧電ｅ定数マトリクスが計算される。
【００５０】
なお、上記形態では、ＰＺＴが正方晶である場合について説明したが、菱面体晶の場合も本発明を適用できる。この場合には、ドメインウォールで区切られているドメイン間の自発分極のなす角度を菱面体晶に対応した角度に変更する必要がある。
【実施例】
【００５１】
先ず、応力、分極およびひずみのうち少なくとも一種の物理量を算出しようとする圧電磁器について、分極した状態での弾性スティフネスマトリクス〔ｃ〕、圧電ｅ定数マトリクス〔ｅ〕、誘電率マトリクス〔ε〕をＥＭＡＳ試験により求め、圧電磁器のドメインの自発ひずみの大きさγｓｐを、Ｘ線回折測定にて求めた。結果を表１に記載した。さらに、圧電磁器のＤ−Ｅヒステリシスカーブ、バタフライカーブを測定した。
【００５２】
次に、本発明の圧電磁器の物理量算出方法を用いた有限要素法プログラムにより、仮想的なＤ−Ｅヒステリシスカーブ、バタフライカーブを求め、図６に示した実測したＤ−Ｅヒステリシスカーブ、バタフライカーブについてフィッティングを行った。フィッティングパラメータを表２に記載した。
【００５３】
この計算は１つの有限要素で行った。また１つの有限要素内では６４９通りの結晶方位についてドメインの動きを計算し、その寄与を合計して有限要素の物理量を計算した。本発明の計算手法では、１つの有限要素で非線形圧電現象を表現することができ、かつ未分極から分極状態への遷移過程についても直線ではなく実験結果を再現できた。
【００５４】
【表１】
2010258188.tif 000008
【００５５】
【表２】
2010258188.tif 000009
【００５６】
この後、多数の要素で計算実行し、応力、分極、ひずみ分布を得た。図１に、応力分布を記載した。なお、表１、２において、例えば、１．２９Ｅ＋１１とは１．２９×１０１１を意味する。
【図１】
2010258188.tif 000010
【図２】
2010258188.tif 000011
【図３】
2010258188.tif 000012
【図４】
2010258188.tif 000013
【図５】
2010258188.tif 000014
【図６】
2010258188.tif 000015