コラム

動的陽解法について

陽解法のすすめ/第2回担当:里 優(2016.2)


前回は陽解法の特徴を説明しました。 今回は、もう少し詳しく動的問題の陽解法(動的陽解法と呼ばれます)について述べてみたいと思います。 なお、離散化などに関する詳しい説明は、「動的陽解法と動的緩和法について」にまとめてありますのでご覧ください。


解くべき運動方程式は、次式です。


式-1

ここに、uは変形、fiは内力、feは外力、ρは密度、cは減衰定数です。


時間微分を中央差分で近似して


式-2

式-3

とすれば、式(1)は、


式-4

となります。


要素分割した領域に対して、内挿関数を重みとして式(4)にガラーキン法を適用すると、要素毎の合算が終わった形での節点値として次式が得られます。


式-5

式-6

ここに、Fsは表面力、Fiは内部応力の等価節点力、Fbは物体力(体積力)です。


前回でも説明しましたが、係数マトリクスの[C]は対角化してあり、上式は連立方程式ではなく節点ごとの並列式となります。 したがって、節点単位で次のようにして変位の更新値を求めることができます。


式-7

このように、動的陽解法では要素単位で節点値を求めて行くため、全体剛性マトリクスを作成することも、連立方程式を解くこともありません。 計算のためのコーディングも極めてシンプルなものとなります(「動的陽解法と動的緩和法について」参照)。


それでは、実際に動的陽解法で問題を解いてみましょう。 なお、今後ご紹介する計算例では、1次の内挿関数で内部変数を定義する2次元の三角形要素を用います。 これは、係数マトリクスの対角化が容易で物理的な意味もわかりやすいと考えているからです。 4角形要素では、図のように要素を4つの三角形要素に分割して節点値などを求め、積算した値を2で割って平均化する方法を用いています。


図-1 四角形要素の取り扱い
図-1 四角形要素の取り扱い

解析モデルを下図に示します。100mの地盤を棒状に切り取ったイメージです。物性値は、表-1にまとめてあります。


図-2 解析モデル
図-2 解析モデル

表-1 解析に用いた物性値一覧
物 性 単位
E ヤング率 160 MPa
ν ポアソン比 0.3
γ 単位体積重量 15 kN/m3
g 重力加速度 9.8 m/s2
Vp P波速度 375 m/s
Vs S波速度 200 m/s

このモデルの底面に、図-3のような強制水平変位を与えました。 なお、各節点の鉛直変位は拘束しています。また、減衰は加えていません。


図-3 底面に加えた強制変位
図-3 底面に加えた強制変位

結果を動画でご覧にいれます。動的陽解法によって変形の経時変化が求められています。また、最上部の節点における水平変位の経時変化を図-4に示します。 Vsが200m/sなので、100mのモデルの最上部に変形が到達するのが0.5秒となりますが、計算結果には多少誤差が認められます。 これは、要素分割や時間刻みの影響と思われます。


図-4 モデル最上部の水平変位
図-4 モデル最上部の水平変位

動画-1 解析結果

動的陽解法で注意しなければいけない点は、時間刻みΔtの制限です。 これは、Courant条件と呼ばれ、次のように表されます。


式-8

式-9

ここに、VpはP波速度、Lminは解析対象を構成する要素の中で最も小さな節点間距離、ρは密度です。


この制限は、Δtの間にP波が次の要素まで到達しないように時間刻みを設定するというもので、これを超えると計算が不安定となります。 このため、小さい要素があったりヤング率が大きな材料を扱う場合には、Δtを小さくとる必要があります。


次回は、動的陽解法における境界条件の処理をご紹介する予定です。