JP2007080095A

JP2007080095A - 異方性を有する球状表面波素子の解析方法

Info

Publication number: JP2007080095A
Application number: JP2005269025A
Authority: JP
Inventors: Ichiji Yamanaka; 一司山中; Yasushi Ikegami; 泰史池上
Original assignee: Tohoku University NUC; Itochu Techno Solutions Corp
Current assignee: Tohoku University NUC; Itochu Techno Solutions Corp
Priority date: 2005-09-15
Filing date: 2005-09-15
Publication date: 2007-03-29

Abstract

【課題】球状弾性表面波素子および球状光表面波素子の実現の為に、異方性単結晶の弾性表面波伝搬特性を計算する有限要素法極座標系波動伝搬コードを動作させる解析装置を提供する。
【解決手段】異方性単結晶の弾性表面波伝搬特性を計算する有限要素法極座標系波動伝搬コードを動作させる解析装置において、極座標系を用いた球殻モデルを作成し、球殻上に表面波伝搬経路を変化させるために設置された障害物を作成し、球表面上に超音波送受信素子であるIDTモデルを作成し、自動メッシュにより極座標を用いた６面体要素を作成し、解析条件に従い有限要素法極座標系波動伝搬コードを実行する。なお本コードは極座標における要素の形状対称性および材料における要素物性値の対称性を利用し計算の高速化を行ったことを特徴としている。また解析結果として球表面および球内部の変位、ひずみまたは応力値のスナップショット、または時刻歴波形を得ることができる。
【選択図】図２

Description

本発明は、球状弾性表面波素子および球状光表面波素子など異方性を有する球体による素子の効率的設計とプロセス開発の為に開発された、異方性単結晶球の表面波伝搬特性を計算する極座標系有限要素波動伝搬コードおよびこれを動作させる解析装置に関するものである。

最近、球状弾性表面波素子など球体の性質を利用した素子が開発され、センサなどの飛躍的高性能化に貢献することが期待されている。従来、この種の素子の解析方法としては、等方性を仮定して成立する解析解を用いた解析方法や、直方体要素を用いる有限要素法が用いられていた。しかし解析解では強い異方性を高精度に扱えない問題があった。また、従来用いられる直方体要素を用いるものは、球の表面を忠実にモデル化できず、精度が低いという問題があった。

解析解を用いる球状弾性表面波素子の解析技術は、例えば非特許文献１に示されているが、図１(a)のように等方性材料を仮定した解を高速かつ精度良く得ることができる優れた技術であるが、弾性表面波素子に良く使われる強い異方性材料ではあまり精度が良くなかった。また、直方体要素を用いる有限要素法は、特許文献１に示されているが、図１(b)のように球表面が階段状になっており、球表面を滑らかに表現することが出来ず、肝心の表面近傍での解析精度が悪かった。

特開２００２−１２２６７５号公報：「地質波動伝搬シミュレーションシステム及びその記録媒体」 Kazushi YAMANAKA,et. Al., Ball SAW Device for Hydrogen Gas Sensor, IEEE ultrasonic symposium-299(2003).

上記のように従来の技術では、解析解では強い異方性を持つ材料での精度が悪く、直方体要素を用いた有限要素法では表面が階段状になり、弾性表面波素子の設計または解析に重要な表面での誤差が大きいという問題があった。

本発明は、極座標を用いて球表面を滑らかにモデル化でき、弾性表面波素子の設計または解析に重要な球表面での誤差を最小限にし、かつ高速に計算できる特徴を持つ球に特化した専用コードを提供することを目的としている。

本発明によれば、異方性を有する球状弾性表面波素子の設計解析のための波動伝搬特性に関する計算を、極座標において要素を定義した有限要素法を用いて行うことを特徴とする方法および装置が得られる。

また本発明は、有限要素法極座標系波動伝搬コードに関し、極座標系の６面体有限要素法による解析機能、等方性弾性材料及び異方性弾性材料を扱える機能、波動伝搬解析機能、波動伝搬解析を陽解法で実施する機能、領域分割法による並列化機能を有することを特徴とする方法および装置を提供する。

また本発明は、ユーザの操作環境を向上させるため、極座標系で球殻モデルおよびメッシュを自動生成する機能、球殻上に表面波伝搬経路を変化させるために設置された障害物をモデル化する機能、球表面上に超音波送受信素子であるIDTをモデル化する機能、変位および応力などユーザが求める出力を設定する機能、全伝搬計算時間、球表面および球内部の変位または応力値のスナップショット出力時間間隔、球表面のすべての要素位置における時系列変位および応力波形出力機能を設定することを特徴とする方法および装置を提供する。

また本発明は、極座標における要素の形状対称性および材料における要素物性値の対称性を利用し、計算の高速化を可能にすることを特徴とする方法および装置を提供する。

本発明によれば、極座標を用いて球表面を滑らかにモデル化できるようにし、かつ極座標における要素の形状対称性および材料における要素物性値の対称性を利用し、計算の高速化を行ったので、現実的な計算時間で強い異方性を持つ球状弾性表面波の高精度計算を実現することが可能になる。

図２は本発明の実施の形態による「異方性を有する球状弾性表面波素子の設計解析装置」の概略構成を示す図である。図２を参照すると、まず、シミュレーションに必要な入力ファイルを作成する（ステップａ）。ここでは、極座標系を用いた球殻モデルの定義、球殻上に表面波伝搬経路を変化させるために設置された障害物の定義、球表面上に超音波送受信素子であるIDTモデルの定義、変位、ひずみ、応力などユーザが求める出力種類の定義、全伝搬計算時間、球表面および球内部の変位、ひずみまたは応力値のスナップショット出力時間間隔、球表面のすべての要素位置における時系列変位、ひずみまたは応力波形出力時間間隔の定義を行う。この入力ファイルから前処理（ステップｂ）によりメッシュファイルおよび中間入力ファイルを作成する。なおメッシュファイルを後処理（ステップｆ）で表示することによりメッシュ形状を確認することができる。これら２つのファイルは有限要素法極座標系波動伝搬コード（ステップｄ）の入力となり、有限要素法極座標系波動伝搬コードからは時刻歴ファイル(表面要素)およびコンターファイル(全要素、表面要素)が出力される（ステップｅ）。後処理（ステップｆ）では、時刻歴波形表示、メッシュファイルおよびコンターファイルの立体図および任意断面、包絡面について可視化をおこなう。

また、本発明で対象としている球殻モデルとしては、実施例に記載した球状弾性表面波素子に限らず、モデル化パラメータを調整することにより、地球のような巨大な球体も対象とできるため、特に、用途を球状弾性表面波素子に限定するものではない。

以下、図２の概略構成図に沿って、本発明の実施の形態について図面を参照しながら説明する。

前処理に関して、極座標系で球殻モデルおよびメッシュを自動生成する機能を有するが、オペレータは、異方性を有する球体のどの部分にメッシュ生成を行うかを、内球半径、外球半径、最小余緯度、最大余緯度、最小経度、最大経度、動径方向の分割数、緯度方向の分割数、経度方向の分割数をパラメータに決定する。なお、この際用いられる極座標系は、先の図３で説明したオイラー角により回転したxyz_e座標系で定義される。また、球殻の物性値は剛性マトリクスの各成分を直接入力することにより決定する。

前処理に関して、球殻上に表面波伝搬経路を変化させるために設置された障害物をモデル化する機能を有するが、球殻上の障害物は、図３に示されるように障害物の外形形状を複数の直線で囲まれた領域とし、動径方向に一様な厚みを持つと定義する。オペレータは、障害物の物性値ＩＤ、外形を定義する節点数(nid)、厚さ(th)を定義し、さらにこれら節点の球殻上の座標値を余緯度、経度(thn_e,phn_e)をパラメータに決定する。

前処理に関して、球表面上に超音波送受信素子であるIDTをモデル化する機能を有するが、オペレータは、図４に示されるようにIDTの幅(widt)、長さ(leng)、厚さ(thick)、設置間隔(pich)、設置数(num)、IDTの始点座標値を余緯度、経度、ねじれ角をパラメータに決定する。また、各振動子に入力する超音波は、それぞれの振動子ごとに振動方向ベクトルの時系列波形を与えることにより実現する。また、ＩＤＴの物性値は剛性マトリクスの各成分を直接入力することにより決定する。

前処理に関して、変位、ひずみ、応力などユーザが求める出力を設定する機能を有するが、オペレータは、変位、ひずみ、応力のどの成分を出力するかの選択、出力する領域を球殻すべての要素とするか球表面のみとするかの選択を設定することができる。

前処理に関して、全伝搬計算時間、球表面および球内部の変位、ひずみまたは応力値のスナップショット出力時間間隔、球表面のすべての要素位置における変位、ひずみまたは応力波形出力設定機能を有するが、オペレータは、全伝搬計算時間、スナップショットまたは時系列波形の出力時間間隔を設定することができる。

有限要素法極座標系波動伝搬コードに関して、極座標系の６面体有限要素法による解析コードであるが、本実施形態では、図５に示すように結晶軸は基準座標系のｚ₀軸方向を向いているとし、６面体要素を生成する領域は灰色で色づけされた部分で、オイラー角により回転したxyz_e座標系での極座標で定義されている。

有限要素法極座標系波動伝搬コードに関して、材料として、等方性弾性材料及び異方性弾性材料を扱えるが、本実施形態では、材料定義には、数式（１）に示される応力−ひずみ関係式で定義される剛性マトリクスＤ（６×６の対称マトリクス）の変数を直接入力することにより実施する。これにより、等方材を含む任意の異方弾性材料を取り扱うことが可能になる。
σ＝Ｄε （１）

有限要素法極座標系波動伝搬コードに関して、陽解法の波動伝搬解析機能を有するが、本実施形態では、有限要素法の定式は次の通りである。

[Ｍ][ｕ2]＋[Ｃ][u1]＋[Ｋ][u]＝[Ｐ] （２）

ここで、[Ｍ]、[Ｃ]、[Ｋ]は、それぞれ質量マトリクス、減衰マトリクス、剛性マトリクスであり、[u2]、[u1]、[u]は、それぞれ加速度ベクトル、速度ベクトル、変位ベクトルである。また、[Ｐ]は、外力ベクトルである。

そこで、数式（１）を時間方向にｄｔの間隔で離散化すると、時刻t＋ｄｔでの変位u_t+dtは、数式（３）で表わされる。

[u_t+dt]=[(ｄｔ)²・[P_t]＋(2[Ｍ]−d t・[Ｃ]−(ｄｔ)²・[Ｋ])・[u_t]＋（d t・[Ｃ]−[Ｍ]）・[u_t-dt]]／[Ｍ] （３）

ここで、集中質量を採用すると、数式（３）の[Ｍ]^-1は、対角マトリクスとなり、陽解法で変位[u_t-dt]を計算することができる。

有限要素法極座標系波動伝搬コードに関して、波動伝搬解析は、領域分割法による並列化コードとするが、本実施形態において、３つの領域で分割した際の領域分割の例を図６に示す。なお、領域間の節点は図７に示すように２重節点として、これらの節点どうしで節点力をＭＰＩによりデータ交換することで並列計算を実現している。ただし、領域分割の例は本実施例では、経度方向に分割したが、データ交換方法に変更なく深さ方向に分割することも可能であり、分割方向を限定するものではない。

有限要素法極座標系波動伝搬コードに関して、極座標における要素を利用し、計算の高速化を可能にする機能を有するが、極座標系で要素を作成すると、要素ごとに方位が異なるので剛性マトリクスも異なるため、立方体要素を用いる場合に較べて、メモリ使用効率や計算効率の面で大幅に不利になる。

本実施形態では、有限要素法の計算で最も時間のかかる節点力計算ループにおいて、有限要素法でメモリを最小化する際に一般的に用いられる各要素で要素剛性マトリクスを計算する手法を採用するが、このままでは計算時間が膨大に必要となるため、新たに図８に示されるアルゴリズムを適用した。これによりメモリ、および計算時間の両方の最小化が可能となった。

すなわち、通常、各要素で要素剛性マトリクスを計算する手法では、形状関数の計算は節点力計算ループ（ステップｂ）の中に入るが、球体の経度方向対称性を利用し、経度０度のみの断面の形状関数を先に計算しメモリに保存する。この際のメモリは断面のみの要素分であるため、メモリ負担は微小である。次に節点力ループ（ステップｂ）に入ると、要素節点変位ベクトルを先に計算した形状関数の対象とする要素まで経度方向に回転する（ステップｃ）。さらに異方性を示す材料マトリクスも同じく経度方向に回転する。この際、材料マトリクスの対称性が大きい、すなわち横方向等方材や完全等方材であれば、回転に要する自由度を減らせるまたは無くすことができるため、高速化を図る事が可能である。次に要素剛性マトリクスを計算する（ステップe）。この計算は、回転後の材料マトリクスとメモリに保存された形状関数の単純なかけ算となるので、計算負荷は通常の積分計算を伴うものより大幅に削減することが可能になる。次に要素節点力の計算を行う（ステップｆ）。その際得られた節点力は経度０度での値であるため、本来の経度Φn度に回転する（ステップｇ）。最後に全体節点力ベクトルへ回転後の要素節点力ベクトルを組み込み（ステップｈ）、次の要素の計算に移る。

後処理に関して、図８に示されるようにメッシュ可視化機能、時系列波形出力機能、３次元可視化機能、任意断面可視化機能、包絡面可視化機能を有する。メッシュ可視化機能では、有限要素に分割された計算対象領域、障害物、およびIDTの形状を画面出力する。時系列波形出力機能では、外球面上の任意の観測点における場の時間発展の様子を調べる。コントロールパネル上で観測点(余緯度と経度)を入力、変位、ひずみ、または応力などのユーザが求める出力を選択するとグラフが画面出力される。３次元可視化機能では、外球面上を伝搬する波動パターンを立体的に可視化する。変位、ひずみ、または応力などのユーザが求める出力を選択すると立体的な波動伝搬パターンの動画を画面出力する。任意断面可視化機能では、球を横切る任意切断面内の波動パターンを可視化する。変位、ひずみ、または応力などのユーザが求める出力を選択し、切断面の設定と切断面の法線方向の平行移動する距離を設定すると切断面が設定され、この断面内の波動伝搬パターンの動画を画面出力する。包絡面可視化機能では、外球面上の波動パターンの包絡面を立体的に可視化する。変位、ひずみ、または応力などのユーザが求める出力を選択すると選択された値の外球面上の包絡面を立体的に画面出力する。

本発明に係る「異方性を有する球状表面波素子の設計解析計算方法および装置」は、燃料電池用の水素センサをはじめに、生体分子センシングの飛躍的な高感度化を実現する球状弾性表面波センサの開発等の様々な分野に適用できる。

また本発明によって製造された極座標系有限要素法コードは、地球を周回する地震波伝搬などにも応用できると考えられ、直径数ミリの素子から直径1万2000Kmの地球まで幅広いスケールの波動伝搬解析に応用することができる。

本実施形態による従来の技術の概念図である。本実施形態による「異方性を有する球状弾性表面波素子の設計解析装置」の概略構成を示す図である。本実施形態による表面波伝搬経路を変化させるために設置された球表面上の障害物モデル化の概念図である。本実施形態による超音波送受信用に設置された球表面上のIDTモデル化の概念図である。本実施形態による極座標系の６面体要素のメッシュ生成領域およびメッシュを示す概念図である。本実施形態による３つの領域で分割した際の領域分割の例を示す概念図である。本実施形態による領域間の節点で節点力をＭＰＩによりデータ交換する概念図である。本実施形態による計算の高速化に関するプログラムフローチャートである。本実施形態によるメッシュ可視化機能、時系列波形出力機能、３次元可視化機能、任意断面可視化機能、包絡面可視化機能に関する概念図である。時系列波形出力と３次元可視化と任意断面可視化と包絡面可視化とを示す図である。

Claims

異方性を有する球体の設計解析のための波動伝搬特性に関する計算を、極座標において要素を定義した有限要素法を用いて行うことを特徴とする方法および装置。
請求項1の方法および装置において、有限要素法極座標系波動伝搬コードに関し、極座標系の６面体有限要素法による解析機能、等方性弾性材料及び異方性弾性材料を扱える機能、陽解法による波動伝搬解析機能、領域分割法による並列化機能を有することを特徴とする方法および装置。
請求項1に記載の方法および装置において、ユーザの操作環境を向上させるため、極座標系で球殻モデルおよびメッシュを自動生成する機能、球殻上に表面波伝搬経路を変化させるために設置された障害物をモデル化する機能、球殻上に超音波送受信素子であるIDTをモデル化する機能、変位および応力などユーザが求める出力を設定する機能、全伝搬計算時間、球表面および球内部の変位または応力値のスナップショット出力時間間隔、球表面のすべての要素位置における時系列変位および応力波形出力機能を設定することが可能であることを特徴とする方法および装置。
請求項1に記載の方法および装置において、極座標における要素の形状対称性および材料における要素物性値の対称性を利用し、メッシュ生成および計算の高速化を可能にすることを特徴とする方法。
請求項１に記載の方法および装置において、球体の経度方向対称性を利用し、まず経度０度の断面の形状関数を計算してメモリに保存し、次に所定経度にある要素節点変位ベクトルを先に計算した形状関数の対象とする要素まで経度方向に回転し、次に異方性を示す材料マトリクスも同じ経度方向に回転し、次に回転後の材料マトリクスとメモリに保存された形状関数のかけ算により、要素剛性マトリクスを計算し、次に経度０度での要素節点力ベクトルの計算を行い、次にこの要素節点力ベクトルを前記所定経度度にもどすように回転し、最後に全体節点力ベクトルへ組み込むことを特徴とする方法。
請求項1、２および３に記載の方法および装置、請求項４に記載の方法を製造工程の一部に含むことを特徴とする装置の製造方法。
請求項６に記載された製造方法で製造されたことを特徴とする方法および装置
請求項１から６までにおいて、波動が弾性波であることを特徴とする方法および装置
請求項１から６までにおいて、波動が電磁波であることを特徴とする方法および装置