時間積分の安定性について (1)安定限界

Altair_Ichikawa
Altair_Ichikawa New Altair Community Member
edited July 2018 in 質問と回答 (Q&A)

衝撃問題等を陽解法コードで解く際、これ以上Δtを大きくすると解が発散するという限界、いわゆる安定限界がある事が知られています。

これは、Courantの安定条件と呼ばれ、次式で表されます。

 eq1-1.jpg.f5a3b03e20bda06daa5a03be673315be.jpg                                     

ここに、は要素の代表長さ、cは物質の弾性波速度で、eq1-2.jpg.d264d27d3f5c512c7b5e5d5ccc9d9e36.jpgE:ヤング率、ρ:密度です。これが全ての要素に対して成り立つ必要がありますので、もし、物体の材料が均一ならば、1時間ステップΔtで弾性波が進む距離が最小の要素の長さよりも小さくなるようにする必要がある、というのがこの式の示す所となります。

ところで、数式の解が発散するかどうか、というのは物理の問題ではなく、あくまでも数学の問題です。言ってみれば、Courantの安定条件は、数学的な安定条件の物理的な解釈からなる近似式と言う事ができます。

それでは、数学的な見地から見た場合、この安定条件は、どのように導かれるのかを以下に見ていきたいと思います。

衝撃の問題であれ、振動の問題であれ、その支配方程式は、運動方程式で表す事ができます。

 eq1-3.jpg.6fb1b9e2bdc9ed90a16eb61ec16a6d9e.jpg 

n自由度の運動方程式は、振動のモード解析で知られているように、

 eq1-4.jpg.25ace20e8e9dd9aead05483267d0a32e.jpg

といった、n個のモード座標系での運動方程式の解の重ね合わせで表現できます。このどれか一つでも発散すれば、解は発散することになりますが、その中で一番条件的に厳しいのは、その最大固有振動数の場合ですので、このケースで解が発散しない事を示せれば、安定限界が求められる事になります。

一般に運動方程式を時間方向に解いていく際、その式は、外力等の既知の項を無視すれば、

  eq1-5.jpg.649d299ecb8d8b164c58c3022591b613.jpg

という様な漸化式で表されます、mステップ先まで計算すれば、結局、

 eq1-6.jpg.9c652462cd56feeca4a29ece75ddf00b.jpg 

といった計算を行っている事に相当します。この係数行列Aが積分作用素行列と呼ばれます。この行列の固有モードマトリックスをP、その固有値をλと置くと、Aが2×2の行列であれば、

  eq1-7.jpg.ab4aaa857775f0aa821860cb018c85c2.jpg

の関係が成り立ちますので、結局、

  eq1-8.jpg.8b4ef0c531b28ee26b379b97ef6d4dc3.jpg

が得られ、mがいくら大きくなっても解が発散しないためには、作用素行列の固有値の絶対値が1を超えない事が必要である事が解ります。即ち、

   eq1-9.jpg.9c36aa55ea2776334c1e4ac442076891.jpg

です(ρはスペクトル半径と呼ばれます)。積分作用素行列Aの具体的な形は、用いた時間積分公式によって決まります。陽解法コードでは中央差分法が用いられますので、それによって安定限界が決まり、更にその近似としてCourantの安定条件が得られていた事になります。

 

   

 

Tagged: