FrontISTR Commons 2024年7月20日目次
1本ソフトウェアは文部科学省次世代IT基盤構築のための研究開発「イノベーション基盤シミュレーションソフトウェアの研究開発」プロジェクトによる成果をシーズとして，継続的に開発されている並列有限変素解析プログラムです。本ソフトウェアを無償または営利目的でご使用になる場合、「MITライセンス」をご了承頂くことが前提となります。

1.1 マニュアルリスト
•イントロダクション•インストールマニュアル•理論マニュアル●解析マニュアル•チュートリアル• FAQ本マニュアルでは、大規模並列FEM非線形構造解析プログラムFrontISTRで用いられる、有限要素法（Finite Element Method）による解析手法について説明します。固休の応力解析手法については、まず微小変形線形弾性静解析手法について示し、引き続き大変形問題を扱う際に必要となる幾何学的非線形解析手法、彈塑性解析手法について示します。さらにFEMによる応力解析の結果を利用して得られる破壊力学パラメータを評価する方法についてまとめたものを示します。次に、固有値解析および熱伝導解析手法について示します。
1.2 本マニュアルの記載内容
• PDF •静的解析－微小変形線形弾性静解析－非線形静解析手法•動的解析—動的解析手法–固有値解析－周波数応答解析• 熱伝導解析•参考文献
1.3 微小変形線形弾性静解析
ここでは微小変形理論に基づく弾性静解析についての定式化を示す。応力・ひずみ関係として線形弾性を仮定している。 1.3.1 基礎方程式固体力学の平衡方程式、力学的境界条件、幾何学的境界条件（基本境界条件）は次式で与えられる（図2.1.1参照)。
ここで、♂は応力、♂は表面力、S は物体力であり、S は力学的境界、S は幾何学的境界を表す。図2.1.1固体力学における境界値問題(微小変形問題)微小変形問題におけるひずみ・変位関係式は次式で与えられる。
また、線形弾性体での応力・ひずみ関係式（構成式）は次式で与えられる。
ここでCは4階の弾性テンソルである。 1.3.2 仮想仕事の原理基礎方程式(1)、式(2)、式(3)と等価である、微小変形線形弾性問題についての仮想仕事の原理は次式のように表される。
さらに構成式(5)を考慮して式(6)は次式のように表される。
式(8)において、εはひずみテンソル、Cは4階の弾性テンソルである。ここで、応力テンソルσとひずみテンソルεを、それぞれベクトル形式でタ、êと表すと、構成式(5)は次式のように表される。
ここで、Dは弾性マトリクスである。ベクトル形式で表された応力δ 、εおよび式(9)を考慮して、式(8)は次式のように表わされる。
式(10)および式(7)が、本開発コードにおいて離散化される仮想仕事の原理である。 1.3.3 定式化仮想仕事の原理式(10)を有限要素ごとに離散化して次式を得る。
要素ごとに、要素を構成する節点の変位を用いて変位場を次式のように内挿する。
(5)このときひずみは、式(4)を用いて次式のように与えられる。
式(12)、式(13)を式(11)に代入して、次式を得る。
式(14)は次式のようにまとめることができる。
ここで、式(16)、式(17)で定義されるマトリクスおよびベクトルの成分は、有限要素ごとに計算し、蛋ねあわせることができる。
式(15)が、任意の仮想変位δUについて成立することにより次式を得る。
一方、変位境界条件式(3)は次式のように表される。
式(18)を拘束条件式(19)のもとで解くことにより、節点変位Uを決定することができる。
1.4 非線形静解析手法
前述したように微小変形問題の解析においては、平衡方程式などの基礎方程式と等価な仮想仕事の原理を用いて、この式を有限変茶により離散化することによって有限変茶解析を行うことができる。構造物の大変形を扱う有限変形問題の解析においても基本的には仮想仕事の原理が用いられる点は同様である。しかしながら、有限変形問題においては、たとえ材料の線形性を仮定しても仮想仕事の原理式は変位に関して非線形な方程式になる。非線形式を解くためには通常、反復法による繰り返し計算が用いられる。その反復計算においては、ある小さな荷重増分に対して区分的に行なわれ、それを積み重ねて最終的な変形状態へと至る増分解析手法が用いられる。微小委形問題を仮定した場合、ひずみや応力を定義するための配置は、変形前と変形後とで5とくに区別を行なっていなかった。すなわち、微小変形を仮定している場合には基礎方程式を記述する配置は変形前であっても変形後であっても問題にはならなかった。しかしながら、有限変形問題において増分解析を実施する場合、参照配置として最初の状態を参照するか、増分の開始点を参照するかの選択が可能である。前者をtotal Lagrange法、後者をupdated Lagrange法と呼ぶ。詳細については章末参考文献などを参照されたい。本開発コードでは、total Lagrange法およびupdated Lagrange法の双方を採用している。 1.4.1 幾何学的非線形解析手法
1.4.1.1 仮想仕事式の増分分解
時刻tまでの状態が既知であり、時刻t  ＝ t + △tの状態を未知とする増分解析を想定する（図2.2.1参照）。静的境界値問題の平衡方程式、力学的境界条件、幾何学的境界条件（基本境界条件）は次の通りである。
ただし‘‘σ, ‘‘Ђ, ‘’π, ‘‘ণ, ‘‘ণは、それぞれ時刻’におけるCauchy応力(真応力)、物体力、物体表面での外向き単位法線ベクトル、既定された表面力、既定された変位である。これらの式は、時刻t’での配置* v, *´s,, *´s,,に対して記述されるものである。図2  増分解析の概念図2.2.1  増分解析の概念1.4.2 仮想仕事の原理式(20)の平衡方程式と式(21)の力学的境界条件と等価な仮想仕事の原理は次式で与えられる。
ここで、ť A zはAlmansiひずみテンソルの線形部分であり、具体的には次式で表される。
式(23)を幾何学的境界条件、ひずみ変位関係式、応力ひずみ関係式とともに解けばよいのであるが、式(23)は時刻t’の配置で記述されており、現段階で時刻t’の配置は未知である。そこで、時刻0の配置Vまたは時刻tでの配置f’vを参照した定式化が行われる。
1.4.2.1  total Lagrange法の定式化
ここでは、開発コードで用いられるtotal Lagrange法に基づく定式化を示す。時刻0の初期配置を基準とする時刻t’での仮想仕事の原理式は、次式で与えられる。
ただしぎS,ぎEは、それぞれ時刻0の初期配置を基準とする時刻ぜでの2nd Piola-Kirchhoff応力テンソル、Green- Lagrangeひずみテンソルを表す。また、お子,おるは、公称表面力ベクトル、初期配置の単位体積あたりに換算した物体力であり、式(20)、式(21)、式(22)と関連させて、次式で与えられる。
時刻tにおけるGreen-Lagrangeひずみテンソルは次式で定義される。
ここで、時刻ぜにおける変位、2nd Piola-Kirchhoff応力ť u, ♂ Sを次式のように増分分解して表す。
このとき、変位増分に関連して、Green-Lagrangeひずみの増分は次式で定義される。式(30)、式(31)、式(32)、式(33)、式(34)、式(35)を、式(25)、式(26)に代入して次式を得る。ここで、ΔSは、ΔE ,と4階テンソルgCと関連づけて次式のように表されると仮定する。式(36)に式(37)を代入し、△uの二次以上の項を有する△S : δ△E NL を省略して次式を得る。式(38)を有限要素により離散化して次式を得る。ここで、♂ K L, ♂K NL, ♂ F, ♂Qは、それぞれ、初期変位マトリクス、初期応力マトリクス、外力ベクトル、内力ベクトルである。したがって、時刻tの状態から、時刻t’の状態を求めるための漸化式は次式で与えられる。 i = 0 (33) 88 Step1 : $ K (0) =$ K L +$ K N L ; $ Q (0) =$ Q; U (0) =$ U
i = i + 1 1.4.2.2    updated Lagrange法の定式化時刻tの現配置を基準とする時刻t 'での仮想仕事の原理式は、次式で与えられる。
ただし
テンソル* S, * Eやベクトル、* 、 *が時刻tの現配置を基準としているが、Green-Lagrangeひずみについては初期変位(時刻tまでの変位) uを含まず
ただし
の形になる。一方
であるから、これを式(40)、式(41)と式(44)に代入し整理すると解くべき方程式が次のように与えられる。 99

1.4.3.1 超弾性材料
ここで、Διδは、ΔιΕιと4階テンソルεCと関連づけて次式のように表されると仮定する。これを式(48)に代入し、次式を得る。式(50)を有限要素により離散化して次式を得る。ここで、↓K L, ↓K NL, ↓ F, ↓Qは、それぞれ、初期変位マトリクス、初期応力マトリクス、外力ベクトル、内力ベクトルである。したがって、時刻tの状態から、時刻t’の状態を求めるための漸化式は次式で与えられる。 i = 0 i = i + 1 1.4.3 材料非線形解析手法本開発コードでは、等方性超弾性および弾塑性二種類の非線形材料を解析することができる。解析で対象とする材料は弾塑性材である場合では、updated Lagrange法を適用し、超弾性材である場合では、total Lagrange法を適用している。また、反復解析手法にはNewton-Raphson法を適用している。以下にこれらの材料構成式の概要を示す。等方性超弾性材料における弾性ポテンシャルエネルギーは、応力の作用していない初期状態からの等方性を持った応答から得られるものであり、右Cauchy-Green変形テンソルCの主不変量( I 1 , I 2 , I 3 )、または体積変化を除いた変形テンソルの主不変量( I , I , I , I 3 )の関数,つまり、W ＝ W ( I , I 2 , I 3 )あるいはW ＝ W ( I , I 2 , I 3 )として表すことができる。超弾性材の構成式は2nd Piola-Kirchhoff応力とGreen-Lagrangeひずみの関係で定義され、その変形解析はtotal Lagrange法を適用する。以下では本開発コードに含まれた超弾性モデルの弾性ポテンシャルエネルギー Wを列挙する。弾性ポテンシャルエネルギー Wがわかれば、以下のように2nd Piola-Kirchhoff応力および応力－ひずみ関係を計算できる
1.4.3.1.1   (1) Neo Hookean超弾性モデルNeo-Hookean超弾性モデルは等方性を持つ線形則(Hooke則)を大変形問題へ対応できるように拡張したものである。その弾性ポテンシャルは以下のとおりである。
ここで、C 10とD 1は材料定数である。 1.4.3.1.2   (2) Mooney Rivlin超弾性モデル
ここで、( 710, C 01とD 1は材料定数である。 1.4.3.1.3    (3) Arruda Boyce超弾性モデル
ここで、μ, λ mとDは材料定数である。
1.4.3.2 弾塑性材料
本開発コードでは、関連流れ則に準じる弾塑性構成式を適用している。また、その構成式はKirchhoff応力のJaumman速度と変形速度テンソルの関係を表し、その変形解析はupdated Lagrange法を適用する。 1.4.3.2.1   (1)弾塑性構成式弾塑性体の降伏条件が次のように与えられるものとする。初期の降伏条件
後続の降伏条件
ここで、 •  F :降伏関数• σ y :初期降伏応力• σ y :後続の降伏応力● σ :応力テンソル• e :微小ひずみテンソル• eP：塑性ひずみテンソル●ぞ:相当塑性ひずみ降伏応力-相当塑性ひずみ関係が、単軸状態での応力-塑性ひずみ関係に一致するものとする。単軸状態での応力-塑性ひずみ関係
ここで、H'は歪硬化係数である相当応力-相当塑性ひずみ関係
後続の降伏関数は一般には温度、塑性ひずみ仕事の関数であるが、ここでは簡単のため相当塑性ひずみでのみの関数であるものとする。塑性変形の進行中はF = 0が満たされ続ける為、次式が成立しなければならない。
式(64)中のFはFの時間導関数を表しており、以後、ある量Aの時間導関数をAで表す。ここで、塑性ポテンシャルΘの存在を仮定し、塑性ひずみ速度を次式で表すものとする。
ここで入は係数である。さらに、塑性ポテンシャルΘが降伏関数Fに等しいものとして、次式の関連流れ則を仮定する。
この式を式(64)に代入し、下式が得られる。
ここで、Dは弾性マトリクスであり、
弾塑性の応力一ひずみ関係式は以下のように書ける。
弾塑性材の降伏関数Fがわかれば、この式からその構成式が得られる。
1.4.3.2.2 (1)降伏関数
以下では本開発コードに含まれた弾塑性降伏関数を列挙する•  Von Mises降伏関数
•   Mohr-Coulomb降伏関数
●   Drucker-Prager降伏関数
ここでは、材料定数αとσ,,は材料の粘着力と摩擦角から以下のように計算する

1.4.3.3 粘弾性材料
本開発コードでは、一般化されたMaxwellモデルを適用している。その構成式は以下のように幅差ひずみeと幅差粘性ひずみqの関数になる。
ここでは
である。また、qは
から求められる。ここではλ,,はリラクゼーションである。また、リラクゼーション係数Gは、以下のProny級数で表す。

1.4.3.4
クリープ材料
応力一定の状況下において時間依存性のある変位は「クリープ」と呼ばれる現象である。前述した粘弾性挙動も一種の線形なクリープ現象と考えることができる。ここでは、いくつかの非線形なクリープの説明を行うこととする。この現象は瞬間的に発生するひずみに追加することで構成式とする方法が一般的に用いられ、ある定荷車が継続している間のひずみをクリープひずみεξ とする。クリープを考慮した構成式は、通常、応力と全クリープひずみの関数として定義されるクリープひずみ速度ες が用いられる。
ここで、瞬間的に発生するひずみが弾性ひずみεξであるとすると、全ひずみはクリープひずみを加えた次式のように表される。
である。前述の塑性材料でも示したように、クリープを示す構成式に対して数値解析上の時間積分の方法を示さなければならい。クリープを考慮したときの構成式は、
ここで、 βπ+θ は、
とする。また、クリープひずみ増分△ες は、非線形方程式を単純化した
とする。 Newton-Raphson法での反復計算では、初期値をσπ+1 ＝σπおよび有限要素法から求められるひずみ増分として、反復解と増分解は次式とする。
ここで、
とする。式(84)と式(85)の解を使って残差Rが0になるまで反復解法を行うとき、応力σπ+1と接線係数
を用いる。式(78)の具体的な式として、本開発コードは、以下のようなNortonモデルを適用している。その構成式は下式のような相当クリープひずみε-rがmises応力qと時間tの関数と表す。
ここでは、A,m,nは材料定数である。 1.4.4 接触解析手法2つの物体が接触すると、接触面を介して接触力むが伝達される。仮想仕事の原理式（2.2.4）を以下のように書きかえる。
ここで、S」は接触面積、u(1)とu(2)はそれぞれ接触物体1と接触物体2の変位を表している。接触解析では、接触する可能性のある面を対にして指定する。この面の対の片方をマスター面、もう片方をスレープ面とする。このマスタースレープ解析手法では、接触拘束条件を以下のように仮定する。 1.スレープ節点は,マスター面を貫通しない。 2.接触があった時、スレープ節点は接触位置とし、この接触点を通じマスター面とスレープ面が互いに接触力，摩擦力を伝達する。式(89)の最後の項を有限要素により離散化して次式を得る
ここでは、K ,とF ,はそれぞれ接触剛性マトリクスおよび接触力を表す。この式を式(39)あるいは(51)に代入すると、接触拘束を考慮したtotal Lagrange法およびupdated Lagrange法の有限要素法定式は以下のようになる。
本開発ソフトは変形体同士問の接触変形解析ができ、ユーザーから以下の解析機能を選択できる。 •微小すべり接触問題：この解析では接触点の位置変化がないと仮定している。 •有限すべり接触問題：この解析は、変形と伴い接触位置変化のある場合に対応している。 16 •摩擦なし接触問題•摩擦あり接触問題：この解析はCoulomb摩擦則に対応している。ただし、微小変形線形弾性解析を選択した場合は、微小すべり摩擦なし問題となる。また、現時点では一次ソリッド要素（要素番号341, 351, 361）の接触解析のみ対応している。
1.5 動的解析手法
本節では直接時間積分法を適用した動的問題解析手法について示す。以下に示すように、本聞発コードでは、陰解法及び陽解法による時刻歴応答解析が可能である。 1.5.1 陰解法の定式化について動的問題を対象として、下式に示す運動方程式の解法に直接時間積分法を適用した。
ここでは、Mは質量マトリクス，Cは減衰マトリクス、Qは内力ベクトル，Fは外力ベクトルである。なお、本ソフトは質量の変化を考慮せず、質量マトリクスは非線形において変形によらず一定となる。時間増分Δt内での変位、速度及び加速度の変化は、Newmark-β法を用いて式(94)及び式(95)に示すように近似している。
ここで、γ, βはNewmark-β法のパラメータである.よく知られているように、γ及びβを以下の値にした場合、線形加速度法あるいは台形則に一致する。
式(94)、式(95)を式(93)に代入すると次式が得られる。
特に、線形問題に対してはK Lは線形剛性マトリクスとし、Q(t + Δt) = K L U (t + Δt)となり、この式を上式に代えすると次式が得られる。
なお、幾何学的境界条件として加速度が指定されている箇所では、式(94)から次式の変位を得る。
同様に、速度が指定されている箇所では、式(98)から次式の変位を得る。
ここで、u is (t + Δt)は時刻t + Δtにおける節点変位、u is (t + Δt)は時刻t + Δtの節点速度、u s (t + Δt)は時刻t + Δtの節点加速度、iは節点自由度番号、sは節点番号である。また、質量項及び減衰項の取り扱いは次のとおりとした。
1.5.1.1質量項の取り扱い
質量マトリックスについては原則として集中質量マトリックスとして扱っている。
1.5.1.2 減衰項の取り扱い
減衰項については式(100)で表されるRayleigh減衰として扱っている。
ここで、R m 、R kはRayleigh減衰のパラメータである。 1.5.2   陽解法の定式化について陽解法では下式に示す時刻tにおける運動方程式を基にする。
ここでは、時刻t + △t及び時刻t − △tにおける変位を時刻§t$におけるTaylor展開により表し、△tに関する2次項までとると、次のようになる。

熱伝導解析
1.6
式(95)、式(96)の差及び和から次式が得られる。式(104)、式(105)を式(101)に代入すると次式が得られる。特に、線形問題に対してはQ(t) = K L U(t)となり、上式は以下になるここで、質量マトリックスMを集中質量マトリックス、減衰マトリックスを比例減衰マトリックスC ＝ R m Mのようにおくと、式(107)は連立方程式の求解操作を不要とする。従って、式(107)からU(t + Δt)は次式により求めることができる。本開発コードで用いられる有限要素法（Finite ElementMethod）による固体についての熱伝導解析手法を示す。 1.6.1基礎方程式連続体中での熱伝導方程式は以下のようになる。ここで、ρ = ρ(x)は質量(密度)、c = c(x,T)は比熱、T = T(x,t)は温度、k = k(x,T)は熱伝導率、Q = Q(x,T,t)は発熱量、zは位置、Zは温度、tは時間を表す。考慮している領域をS、その周囲をΓとする。Γ上では、Dirichet型かNeumann型のいずれかの境界条件が、いたるところで与えられるものと仮定すると境界条件は以下のようになる。
ただし、T 1 ,  qは関数形が既知とする。qは境界からの流出熱流束である。本プログラムでは、3種類の熱流束が考慮できる。
ここで、q,は分布熱流束、q,は対流熱伝達による熱流束、q,は輻射熱伝達による熱流束である。ただし、Tc = Tc(x,t)は対流熱伝達率雰囲気温度、hc = hc(x,t)は対流熱伝達係数、Tr = Tr(x,t)は輻射熱伝達率穿囲気温度、hr = εσF = hr( x, t )は輻射熱伝達係数、εは輻射率, σはStefanBoltzmann定数、Fは形態係数である. 1.6.2 離散化方程式(109)をGalerkin法によって離散化すると、
ただし、
20
方程式(116)は非線形かつ非定常の方程式である。いま、時間に関して後退オイラー法により離散化して、時刻t = toにおける温度が既知のとき時刻t = to+∆tでの温度を次式を用いて計算することにする。
ここでの式(121)を近似的にみたす温度ベクトルT(i)そのために、まず、温度ベクトルを次のようにあらわす。
熱伝導マトリクスと温度ベクトルとの積、質量マトリクスなどを次式のように近似的にあらわす。
式(122)、式(123)、式(124)を式(121)に代入して二次以上の項を省略すると次式を得る。
さらに左辺の係数マトリクスを次式を用いて近似評価する。
ここでK ( i )       は接線剛性マトリクスである。結局次式を用いて反復計算を行うことによって時刻t = t o + Δtでの温度を計算することができる。
特に定常解析においては次式を用いて反復計算を行う。
非定常解析において時間増分$ t$の選び方は、時間に関する離散化に陰解法を採用しているので、一般にその大きさの制約を受けない。ただし時間増分Δtが大きすぎると、反復計算における収束回数は増加する。そこで本プログラムは、反復計算過程における残差ベクトルの大きさをつねにモニターし、反復計算の収束がおそすぎれば時間増分Δtを減少させ、反復計算回数が少なくなると時間増分Δtを増加される自動増分機能を備えている。
1.7 固有値解析
1.7.1 一般化固有値問題連続休の自由振動解析を行う場合、空間的離散化を行い、図2.3.1に示すような集中質点による多自由度系でモデル化さられる。減衰のない自由振動問題の場合、支配方程式（運動方程式）は以下のとおりである。
ただし、uは一般化変位ベクトル、Mは質量マトリックス、Kは剛性マトリックスである。ところで、固有角振動数をωとし、a、b、cを任意定数、xをベクトルとして、関数
を定義する。ここで、この式と、この2階の微分、すなわち、
を式(129)に代入すれば、 Mu + Kü = (asinωt + bcosωt)(−ω²M + Kx) = (−λM + Kx) = 0
となる。すなわち、
を得る。つまり、方程式(133)を満たす係数入＝ω ²およびベクトルェを見つけられれば、関数u(t)は、方程式(129)の解となっている。係数λを固有値、ベクトルェを固有ベクトルと呼び、これらを式(129)から求める問題を一般化固有値問題と呼ぶ。図2.3.1  減衰のない自由振動の多自由度系の例図３　減衰のない自由振動の多自由度系の例
1.7.2 問題設定
式(133)は任意の次数に拡張でき、多くの場面で登場する。多くの物理問題を扱う上では行列はエルミート（対称）であることが多い。即ち、複素行列においては、転置行列が共役複素数になっており、実行列においては対称行列である。つまり、行列Kのij成分をk ijとした時、kの共役複素数をkとおけば、
の関係にある。このマニュアル内では、行列は対称で正定値を仮定する。正定値とは固有値がすべて正、言い換えれば下記の式(135)を常に満足する行列である。
1.7.3　シフト付逆反復法有限要素法による構造解析では、実用上、全ての固有値は必要とせず、高々数個の低次の固有値で十分な場合が多い。ところで、HEC-MWでは大規模な問題を扱うことを想定しており、行列はサイズが大きく非常に疎（署要素が多い）である。したがって、この事を念頭に低次のモードの固有値を効率よく求めることが重要である。因有値の下限をσとした時、式(133)を次式のように変形する（数学的には等価な式である)。
この時、計算に当たっては次のような好都合な性質がある。 1.モードが反転している。 2. ρ周辺の固有値が最大化されている。実際の計算では最大固有値が最初に求まることが多い。そのため主要な収束計算を式(133)よりむしろ式(136)に適用し、ρ周辺の固有値から求めることを狙うものとする。この手法は、シフト付逆反復と呼ばれている。 1.7.4 固有値解法のための算法古典的な方法ではJacobi法がよく知られている。この方法は、行列サイズが小さく密行列である時、有効である。しかしながら、HEC-MWで扱う行列は大規模で疎であるため、この方法は採用せずランチョス(Lanczos)反復解法を採用している。
1.7.5ランチョス法
1950年台にC. Lanczosにより提案されたこの手法は、行列を3重対角化する計算算法であり、下記のような特徴を有している。 •反復収束解法であり、行列を疎のまま計算を進めることができる。 •算法は行列、ベクトル積が中心となっており並列化に適している。 •有限要素メッシュに伴う幾何学的領域分割法に適している。 •求める固有値の個数やモード範囲を限定して効率よい計算を行える。ランチョス法は、初期ベクトルからスタートして順次直交ベクトルを作成し部分空間の基底を求める計算を行うものである。この方法は、別の反復解法であるサプスペース法より高速であると言われ、有限要素法プログラムにて広く使われている。しかしこの手法では、計算機の誤差の影響を受けやすく、ベクトルの直交性が損なわれ、途中で破綻する恐れを避けられない。そのため誤差に対する対策は不可欠である。 1.7.6ランチョス法が持つ幾何学的意味式(136)を次のように変数変換することにより
問題を書き直すと
を得る。適当なペクトルgoに対して行列Aによる一次変換を行う（図2.3.2参照）。 24
図2.3.2  行列Aによるqの一次変換
変換されたベクトルは、元のベクトルとつくる空間の中で直交化される。すなわち、図2.3.2のようないわゆるグラム・シュミットの直交化を行う。そうして得られたベクトルをr1としてそれを正規化（長さ1に）してq1を得る（図2.3.3)。同様な算法によりq1からq2を得る。このときq2はq1, q0両方に直交している（図2.3.4）。同様の計算を続けると互いに直交するベクトルが最大行列の次数まで求まる。図2.3.3    qoに直交なベクトルq1図2.3.4    q1とq0に直交なベクトルq2特にランチョス法の算法はベクトル列Aq 0 , Aq 1 , Aq 2言い換えてAqo, A 2 qo, A 3 qo, ,A n qoに対するグラム・シュミットの直交化である。このベクトル列をKrylov列と呼び、それがつくる空間をKrylov部分空25間とよぶ。この空間においてグラム・シュミットの直交化を行うと、直近の２つのベクトルを用いることによりベクトルが求まる。これをランチョスの原理と呼ぶ。 1.7.7 三重対角化上記繰り返しの中でi+1番目の計算は
と表せる。ただし、
である。これを行列表記すると
となる。ここで、
である。すなわち、式(141)で得られる3重対角行列について固有値計算を行うことにより固有値が得られる。 1.8周波数応答解析1.8.1 定式化減衰を考慮しない場合の周波数応答解析の運動方程式は式のようになる。
これを固有モードごとに展開すると、
となる。これを式(143)に代入すると、を得る。この因有振動数が実数になることを以下のように証明する。ω 3 ＝ λ 3とおき、式(145)の複素共役をとると式(146)を得る。これにU jをかけると、を得る。式(147)より、となる。ここで、質量マトリクスは正定値対称なのでゼロベクトルでない固有ベクトルについてが成り立つ。ゆえに、となり、ω 3 ＝入3は実数となる。ここで、2つの造うモードについて考える。これより、が得られ、固有値が異なる場合には、となる。すなわち異なる固有モードは質量マトリクスについて直交する。同じモードについては質量マトリクスについて規格化する(式(154))ことによって取り扱いが容易になる利点がある。次に、減衰を考慮した場合の周波数応答解析についての定式化を示す。対象とする運動方程式を式(155)に示す。
ここでの減衰項はRayleigh型を想定して(156)のようにして表せるものとする。
因有値解析で得られた因有ベクトルにより変位ベクトルは時刻もにおいて式(157)のように展開できる。
このとき、外力項が調和振動子式の
場合においてb,(t)を決定する。運動方程式(155)は強制振動の形になるので
が成り立つ。b,(t)の展開係数の実部および虚部を求めると式(160) iおよび式(161)の形になる。
となる。
1.9