構造物の振動に関する特性を理解する上で、その構造物の固有振動周波数あるいは固有振動モードを知ることは極めて有効である。有限要素法による解析を行う場合、剛性行列をK、質量行列をMとしてKx＝Mxで与えられる固有値問題を解くことによってその構造物の固有振動周波数あるいは固有振動モードを決めることができる。

　近年、計算機の能力の向上により大規模な数値解析が実用的に実施されるようになってきた。さらに、並列計算も実用化されるに至り、超大規模(〜1000万自由度)な数値解析を実施することも可能となってきた。さらに、メッセージ・パッシング・ライブラリーの充実により、並列計算は専用の並列計算機のみならずワークステーションやパソコンをネットワークにより接続した仮想並列計算機上でも実施することが可能となっている。本研究はこれらの背景をふまえ、並列計算により大規模な自由度をもつ構造物の固有振動周波数あるいは固有振動モードを有限要素法により解析することを主目的とする。

　本研究では2つのアルゴリズムを取り扱った。ひとつは、ランチョス法による求解である。ランチョス法は規模の大きな固有値問題でも収束が早く、高い精度の解が安定に得られるアルゴリズムである。このアルゴリズムを並列化することにより、単純形状の構造物の固有振動周波数及び振動モードを高い精度で、最大100万自由度程度の規模の解析を実施することを試みた。さらに、本システムの実用例として水晶振動子の設計を試みた。

　ふたつめは部分構造合成法による求解である。部分構造合成法は全体モデルの固有値問題を直接解く必要がないので、少ない計算量あるいはメモリ量で大規模な自由度をもつモデルの振動解析が実施できるという特徴をもつ。従来のアルゴリズムをさらに改良することによって、一般形状構造構造物の固有振動周波数及び振動モードを適度な精度で、最大1000万自由度規模の解析を実施することを試みた。また、本システムによる一般形状構造物の解析例として工具(金槌)の試みた。

　ランチョス法を大規模な問題に適用する場合の問題点は行列を三角分解する際の計算量と、三角行列をメモリ上に保存するためのメモリ量である。これらの問題を解決するために、行列の三角分解の並列化アルゴリズムとメモリ節約型連立一次方程式求解法を提案した。提案した三角分解の並列化アルゴリズムを用いることにより、超並列計算機のみならず、ワークステーション・ネットワーク上でも80%以上の高い並列効率を実現した。また、提案したメモリ節約型連立一次方程式求解法を用いることにより必要なメモリ量を10%程度に節減することを可能とした。

　図1にこれらの手法を採用した、ランチョス法を用いた求解システムのフローを示す。プロセッサーの中の1個を親プロセッサーとして、残りを子プロセッサーとする。システム行列K及びMの生成は各子プロセッサーで行い、各子プロセッサーはシステム行列の割り当てられた列のみをメモリ上に記憶する。まず、親プロセッサーは初期ベクトルx₁を決定して子プロセッサーに転送する。続いて、子プロセッサーは行列MをB^TBのように三角行列の積に分解する作業を並列に実施しながら、各子プロセッサーで計算されたB^Tの要素に関するx₂＝B^Tx₁を計算する。各子プロセッサーにおける計算結果を親プロセッサーに転送する。親プロセッサーは子プロセッサーの計算結果を足しあわせて、中間ベクトルとして子プロセッサーに転送する。次に、子プロセッサーは行列KをL^TDLのように三角行列の積に分解する作業を並列に実施しながら、前進代入によりx₂＝DLx₃を解く。一回目の反復では、Lの各ブロックの先頭部分をメモリ上に記憶する。各子プロセッサーにおける計算結果を親プロセッサーに転送する。親プロセッサーは子プロセッサーの計算結果を足しあわせて、次の中間ベクトルとして子プロセッサーに転送する。続いて、子プロセッサーは行列Lの再生を行い、後退代入によりx₄＝L^Tx₃を解く。続いて、子プロセッサーは行列Mの三角分解を実施しながら、各子プロセッサーで計算されたBの要素に関するx₅＝Bx₄を計算する。最後の反復では、Bの各ブロックの先頭部分をメモリ上に記憶する。各子プロセッサーにおける計算結果を親プロセッサーに転送する。固有値の決定手順は第一のシステムと同じである。固有ベクトルの計算は、子プロセッサーで行列Bの再生を行い、後退代入により行う。計算結果は親プロセッサーに転送されて、解として出力される。

図1 ランチョス法を用いた解析システムのフロー

　本手法により100万次元、半バンド幅3372の連立一次方程式を32個のプロセッサー(Hitachi SR2201)により約4.2時間で解いた。この問題を解く際に使用したメモリ量は266Mbytesで、本手法を用いない場合の約9%のメモリ量であった。

　図2に本システムを用いた水晶板の厚みすべり振動の変位図を示す。パラメータを変更して計算することによって水晶振動子の設計などに用いることができる。

図2 水晶板の厚みすべり振動の変位

　図3に示すモデルを例に提案する部分構造合成法を説明する。提案する手法では部分構造の振動モードと境界部の振動モードの線形和により全体の振動モードを近似する。

図3 解析モデルと部分構造への分割

　まず、部分構造の振動モードは各部分構造の境界を拘束したモデルの振動モードを用いる。この振動モードをi(i＝1,3,5,7)で表す。ただし、部分構造A、B、C、Dの非境界自由度を添え字の1、3、5、7で表すものとする。境界部の振動モードは各境界について、それに隣接する2つの部分構造を接続した構造により決定する。例えば境界A-Bに関するモードは部分構造AとBを接続し、これ以外の境界を拘束したモデルにより決定する。これより得られるモードを非境界自由度、境界自由度に分けてi^(AB)(i＝1,2,3)と記述する。同様に境界B-Cに関するモードは部分構造BとCを接続したモデルより求めi^(BC)(i＝3,4,5)と記述する。境界C-Dに関するモードは部分構造CとDを接続したモデルより求めi^(CD)(i＝5,6,7)と記述する。全体の振動モードを近似する固有値問題は

　となる。ただし、

　である。

　さらに本アルゴリズムを、図4に示すようにモデルを二方向に分割した場合に適用する方法を示す。図4に示すモデルを一旦、一方向のみに部分構造に分割した状態を考える(図5)。この場合、上記アルゴリズムを適用することが可能で、部分構造の振動モードとして構造A、B、C、Dの振動モード、境界部の振動モードとして構造AとB、BとC、CとDをそれぞれ接続した構造の振動モードを用いることによって全体モデルの振動モードを計算することが可能である。さらに、構造A、B、C、Dの振動モード、構造AとB、BとC、CとDをそれぞれ接続した構造の振動モードを求めるために、それぞれにもう一度上記アルゴリズムを適用する。このように提案したアルゴリズムを2段階に適用することによって、モデルを二方向に分割した場合にも適用可能であると考えられる。同様に3段階に適用することによって、モデルを三方向に分割した場合にも適用可能であると考えられ、したがって、本アルゴリズムはモデル規模が極めて大きい場合や、形状が複雑な場合にも適用可能であると考えられる。

図表図4 モデルを二方向へ分割した様子 / 図5 モデルを一旦一方向に分割した様子

　提案したアルゴリズムを検証するため、30万と800万自由度モデルの解析を行った。まず30万自由度モデルは寸法1000x100x1000mの直方水晶板で、全体のモデルを100x10x100個の8節点6面体要素に分割した。全体モデルの自由度数は336,633である。このモデルを均等に4x4個の部分構造に分割し、拡張アルゴリズムを適用した。同様に、全体モデルの固有値問題をランチョス法により解き、固有振動周波数を求めた。合成に用いたモード数は部分構造振動モード15個、2個の部分構造を接続した構造の振動モードを30個、4個の部分構造を接続した構造の振動モードを55個用いた。

　得られた固有振動周波数の精度を図6に示す。横軸は振動モードの次数、縦軸は精度を表す。解の精度は全体モデルの固有値問題をLanczos法により解いたときの固有振動周波数を基準とした。得られる固有振動周波数の精度は1.0x10^-6〜1.0x10^-4と高次の振動モードまで良好である。

　800万自由度モデルは30万自由度モデルと同じ部分構造を20x20個接続したモデルである。モデルの正確な自由度数は8,283,033である。このモデルでは、全体モデルの固有値問題をLanczos法により解くことができないので、固有ベクトルの精度のみを調べた。結果を図7に示す。固有ベクトルの精度は1MHzから17MHzまでの1MHzごとに最近接の固有振動周波数の対応する固有ベクトルの精度を調べた。横軸は周波数、縦軸は次式で定義される精度である。

図表図6 固有振動周波数の精度 / 図7 固有ベクトルの精度

　固有振動周波数の精度は高く、固有ベクトルの精度はモデル規模の増大に伴って低下しないという点から、本アルゴリズムは大規模かつ任意形状の構造物の振動解析に適用可能でかつ有効であると考えられる。

　図8に本アルゴリズムによって工具(金槌)の振動モードを計算した結果を示す。解析の自由度数は約100万である。図で色がまだらになっているのは部分構造への分割を表している。このように本アルゴリズムは大規模一般形状構造物の解析の適用することが可能である。

図8 部分構造合成法による工具(金槌)モデルの解析

　工業製品を有限要素法などにより数値的に解析し、その特性を予測することは、より安全で経済的で、かつ優れた性能をもつ製品を設計する上で極めて有効である。しかし、実際の解析対象は一般に複雑で、さらに十分な解の精度を確保しようとすると、解析の規模は極めて大きなものとなる。幸いにも近年、プロセッサー単体の速度向上に加え、並列計算機が実現されたり、メッセージパッシングライブラリーや並列言語が整備され、実質的に大規模な計算能力を手軽に利用できるようになっている。

　さて、すべての構造物は固有振動を有しており、構造物自体が動いたり、なんらかの外力が加えられると振動が励振されることは避けられない。この場合、固有振動周波数及び固有振動モードの決定は構造物の共振特性のみならず、強制振動や応答特性を計算する上で必要である。

　このような背景をもとに、本研究は並列計算により構造物の大規模固有値解析を行うための並列化手法及びアルゴリズムを対象としている。

　本研究ではまず第一にランチョス法というアルゴリズムに注目し、その並列化手法を提案している。この手法は、専用の並列計算機上のみならず、ワークステーション・ネットワーク上でも高効率で動作するという特徴をもっている。さらに大規模な問題を少ないメモリで解析するための省メモリ手法を提案している。これらの手法を適用した解析システムを構築し、10万〜100万自由度規模の解析を高効率で行っている。また、実際に水晶振動子の解析に適用し、エネルギー閉じ込め現象の解明、新規構造の設計、寸法設計などを行っており、この分野への貢献も大きい。

　さらに、大規模な解析を実施するために部分構造合成法という手法に注目し、新規アルゴリズムを提案し、最大1000万自由度規模の解析を実施している。また、提案したアルゴリズムは従来の手法より計算量が少なくかつ解の精度が高いことを示している。また、解析事例として、100万自由度規模で一般形状の解析を実施している。

　本論文は6章により構成されている。第1章では工業製品の設計の現状、計算機の発展の状況、構造物の固有値解析の意義を明確にし、本研究の背景と目的を示している。第2章では従来の固有値解析手法及び解析事例を紹介した上で、従来の手法の問題点及び本研究の方向性を示している。第3章ではランチョス法の並列化手法、メモリ節約手法を提案し、また、実際に専用の並列計算機のみならず、ワークステーション・ネットワーク上でも高効率で動作する解析システムを構築している。本システムにより10万から100万自由度規模の解析を行い、計算時間、並列効率、解の収束性などに考察を加えている。第4章では100万から1000万自由度規模の解析を実現するために部分構造合成法を提案し、数値事例によりその有効性を立証している。第5章では解析事例として、ランチョス法を用いた解析システムによる水晶振動子の実際の設計事例を示している。また、提案した部分構造合成法のアルゴリズムによる大規模任意形状構造物の解析例も示している。第6章では本研究の主要な結論を述べ、将来への展開を示している。以上のように、本論文では大規模な自由度をもつ構造物の固有値解析を行うためにランチョス法の並列化手法を提案し、解析システムを構築し、さらに一般形状に拡張するため新規部分構造合成法のアルゴリズムを提案しその有効性を立証している。したがって、本研究は計算力学および計算機援用工学に寄与するところが大きい。また、構築したシステムを実際に水晶振動子の設計に有効に利用しており、この分野への貢献も評価できる。

　よって本論文は博士(工学)の学位請求論文として合格と認められる。