有限要素法の理論

Contents

• 内挿関数とガウス積分
• 仮想仕事の式による要素剛性方程式の標準的な形式の導出
• ２次元弾性骨組構造解析
• ２次元弾性トラス構造解析
• 格子桁構造解析
• ３次元弾性骨組構造解析
• ２次元弾性応力解析
• 軸対称弾性応力解析
• ２次元飽和-不飽和浸透流解析
• ２次元非定常熱伝導解析
• ２次元弾性骨組構造振動解析
• ２次元弾性骨組構造幾何学的非線形解析
• ２次元弾性骨組構造座屈解析
• １次元非定常熱伝導解析

内挿関数とガウス積分

要素内任意点の物理量を節点物理量で表現する

四角形要素を用いた有限要素法では，要素内任意点の物理量を，内挿関数 $[\boldsymbol{N}]$ を用いて，要素を構成する節点の物理量で表す必要がある．ここでは，要素での積分にはガウス積分を用いることとし，$[−1,1]$ の範囲を有する正規化パラメータ座標 $(a,b)$ を用いて，物理量を表現することを考える．

以下に，各種２次元問題における，要素内任意位置での物理量を，節点での物理量で表した事例を示す．

２次元弾性問題

四角形要素内任意点の変位 $u$,$v$ を，節点変位 $(u_i,v_i)$，$(u_j,v_j)$，$(u_k,v_k)$，$(u_l,v_l)$ を用い以下のように仮定する．

要素内任意点の変位を $\{\boldsymbol{u}\}$，節点変位を $\{\boldsymbol{u_{nd}}\}$ とすれば，

２次元浸透流問題

四角形要素内任意点の全水頭 $\phi$ を，節点水頭 $h_i$，$h_j$，$h_k$，$h_l$ を用い以下のように仮定する．

要素内任意点の水頭を $\phi$，節点水頭を $\{\boldsymbol{h}\}$ とすれば，

２次元熱伝導問題

四角形要素内任意点の温度$T$を節点温度 $\phi_i$，$\phi_j$，$\phi_k$，$\phi_l$ を用い以下のように仮定する．

要素内任意点の温度を $T$，節点温度を $\{\boldsymbol{\phi}\}$ とすれば，

内挿関数とその微分

実際の計算においては，内挿関数の数式表現とその微分が必要となるため，これらを求めておくことにする．

以下，$(a,b)$ は，$a=[−1,1]$，$b=[−1,1]$ の範囲を有する正規化パラメータ座標である．

内挿関数の数式表現

内挿関数の微分

ここで，内挿関数 $N_i$ に関する微分は，Jacobi行列 $[J]$ を用いて以下の通り表せる．

上記より，

以下に$[J]$の要素を示す．ここに$x_{i,j,k,l}$および$y_{i,j,k,l}$は要素を構成する節点座標である．

以上により，$[J]$ の各要素が定まり，$[\boldsymbol{N}]$ の微分値も定めることができる．

なお，この段階では $[J]$ および $[\boldsymbol{N}]$ とその微分値は，変数 $a$ および $b$ の関数として定義されていることに注意が必要である．

ガウス積分

積分は以下に示すガウス積分により行う．

面積積分

積分点を 4 点 ($n=2$) とすれば，$a$，$b$ および重み $H$ の値は下表の通りであり，$a$，$b$ を変化させた 4 回の値の総和が積分の近似値となる．また，重みが $H=1$ であることから，計算は更に単純になる．

線積分

熱伝導解析では，熱伝達境界は要素の辺で指定されるため，線積分を行う必要がある．その考え方を以下に示す．

正規化された座標系 a-b において，節点 i-j 間を辺 [1]，節点 j-k 間を辺 [2]，節点 k-l 間を辺 [3]，節点 l-i 間を 辺 [4] とする．例えば，辺 [1] において境界条件が指定された場合，内挿関数 [N ] において b = −1 とし，面積積分の場合の表 の a = s1 および a = s2 の [N] 値の総和を求めればよいことになる．

線積分の場合の積分点座標 s および重みは下表の通りである．

仮想仕事の式による要素剛性方程式の標準的な形式の導出

仮想仕事の式による要素剛性方程式の標準的な形式の導出を行う．

一般に，外力を受ける物体が静的釣り合い状態にある場合，物体内部の応力による仮想仕事と表面力などの外力による仮想仕事は等しい関係にある．この関係を式で表すと以下の通りとなる．

ここで，物体を１つの有限要素とし，要素内任意点のひずみ$\{\boldsymbol{\epsilon}\}$および要素内任意点の変位$\{\boldsymbol{u}\}$を，以下のように節点変位$\{\boldsymbol{u_{nd}}\}$で表せたとする．

ここに，$[\boldsymbol{B}]$はひずみ-変位関係行列，$[\boldsymbol{N}]$は内挿関数行列である．

また以下に示す線形弾性体の応力-ひずみ関係を導入する．

$[\boldsymbol{D_e}]$は応力-ひずみ関係行列，$\{\boldsymbol{\epsilon_0}\}$は温度変化などによる初期歪である．

次に物体力$\{\boldsymbol{F}\}$は，内挿関数と節点での慣性力$\{\boldsymbol{w_{nd}}\}$を用いて以下のように表せるものとする．

上記関係を用いることにより，仮想仕事式の左辺（内力による仮想仕事）は以下のようになる．

また，仮想仕事式の右辺（外力による仮想仕事）は以下のようになる．

以上により，次に示すような，一般的な剛性方程式の形が得られる．

上式は，変位を求めるための式であり，初期歪は荷重項として扱っている．要素内部応力は，前出式

を用いて計算する必要があることに注意する．ここに，$\{\boldsymbol{\epsilon}\}$ は，剛性方程式を解いて求めた変位から求めた要素の歪ベクトルである．

２次元弾性骨組構造解析

要素剛性方程式

要素剛性方程式は，以下に示すとおりである．

要素自由度	要素座標系（x-y）および全体座標系（X−Y）

座標変換行列

全体座標系から要素座標系への座標変換マトリックスは以下に示すとおりである．

全体剛性方程式

系全体で組み立てた全体剛性方程式は，以下のようになる．

ここで，要素節点温度荷重ベクトル $\{\boldsymbol{f_t}\}$ は部材の軸力の補正項なので，全体剛性方程式に組み込む場合，座標変換が必要であるが，節点慣性力ベクトル $\{\boldsymbol{f_b}\}$ は，加速度の方向が全体座標系で定義されているので，全体剛性方程式に組み込む場合，座標変換は不要であることに注意する．

また，全体剛性方程式の解としての変位から部材断面力を算出する場合，軸力については以下の温度変化の補正を行う必要が有ることに注意する．

２次元弾性トラス構造解析

要素剛性方程式

２次元トラス構造解析の要素剛性は，２次元骨組構造解析用要素剛性から，モーメント・回転量に関する自由度を除いたものになる．座標変換行列も同様に，モーメント・回転量に関する項を除いたものになる．

２次元トラス構造解析用要素剛性方程式は，以下に示すとおりである．

座標変換行列は，以下に示すとおりである．

格子桁構造解析

要素剛性方程式

格子桁構造解析の要素剛性と２次元骨組構造解析用要素剛性の関係は，以下の通りで，剛性行列を書き換えることにより，相互に利用可能となる．

格子桁構造解析用の要素剛性方程式は，以下に示すとおりである．

座標変換はX-Y平面内で行われるため，座標変換行列は，２次元骨組構造解析のものと同一となる．

全体座標系	要素座標系と要素自由度

(参考)充実矩形断面のねじり定数 $J$

コンクリートなどの充実矩形断面のねじり定数 $J$ は以下により計算できる．

ここで，次式に示す係数 $\eta$ を導入すれば，係数 $\eta$ と部材寸法 $h/b$ の関係は以下の通りとなる．

上記計算は，以下に示す Python プログラムで行った．

import numpy as np

### Result ###
# h/b  1.0    2.0    3.0    5.0    10.0   20.0
# eta  7.114  4.373  3.798  3.433  3.202  3.098
# J=h*b**3/eta

n=100
b=1.0
#hh=np.arange(1.0,21.0,1.0)
hh=np.array([1.0,2.0,3.0,5.0,10.0,20.0])
for idx,elem in enumerate(hh):
    h=elem
    s=0.0
    for i in range(1,n+1):
        s=s+1/(2*i-1)**5*np.tanh((2*i-1)*np.pi*h/2/b)
    x=1/3*(1-192/np.pi**5*b/h*s)
    print('{0:10.3f} {1:10.3f}'.format(h/b,1/x))

３次元弾性骨組構造解析

要素剛性方程式

節点変位ベクトル

節点力ベクトル

剛性行列

座標変換行列

部材軸（x軸）の方向余弦

• １つの梁要素の節点を $i$ および $j$ とする．
• 節点 $i$ の座標を， $(x_i,y_i,z_i)$, 節点 $j$ の座標を $(x_j,y_j,z_j)$ とする．
• この時，部材軸（x軸）の方向余弦は以下の通り表現できる．

座標変換行列の小行列

$\boldsymbol{[t]}$ は，全体座標系を部材座標系に変換する座標変換行列 (12 x 12) の小行列 (3 x 3) である．

ここに， $\boldsymbol{\{x\}}$ は部材座標， $\boldsymbol{\{X\}}$ は全体座標 である．

(1) 部材長軸（部材座標系x軸）が全体座標系Z軸と平行でない場合

(2) 部材長軸（部材座標系x軸）が全体座標系Z軸と平行な場合

座標変換行列

コードアングルが０の場合の全体座標系と部材座標系の関係

２次元弾性応力解析

有限要素式

節点変位計算式

要素応力算定式

節点変位計算式では，温度変化による項 $\{\boldsymbol{f_t}\}$ を荷重ベクトルに加えて変位を計算しているが，要素応力を求める際は，上式に示すように，節点変位から求めた歪から温度変化による歪を差し引く必要が有ることに注意が必要である．

節点変位ー任意位置ひずみ関係式

節点変位ー任意位置歪算定式

２次元弾性問題における応力-ひずみ関係

一般に３次元問題における等方均質弾性体の応力-ひずみ関係は，弾性係数を$E$，ポアソン比を$\nu$，材料の線膨張係数を$\alpha$，温度変化量を$T$として以下の通りとなる（Hookeの法則）．

ここで$x-y$平面上を考える場合，平面応力状態では，

であるため，

平面ひずみ状態では，

であるため，

となる．これらの関係を再整理すると以下の通りとなる．

２次元問題での応力とひずみの成分

応力・ひずみの符号は，引張を正とし，温度の符号は温度上昇を正とする．

２次元問題での応力-ひずみ関係式

4節点4ガウス積分点パラメトリック要素としての定式化

ひずみ-変位関係行列$[\boldsymbol{B}]$の導入

四角形要素内任意点の変位$u$,$v$を以下のように仮定する．ここに座標($a$,$b$)は[$-1$,$1$]の範囲を有する正規化パラメータ座標，$u_{i,j,k,l}$および$v_{i,j,k,l}$は要素を構成する節点の変位を示す．

要素内任意点の変位を$\{\boldsymbol{u}\}$，節点変位を$\{\boldsymbol{u_{nd}}\}$とすれば，

要素内任意点のひずみ$\{\boldsymbol{\epsilon}\}$は，節点変位$\{\boldsymbol{u_{nd}}\}$を用いて，

要素剛性行列

4節点パラメトリック要素の剛性行列$[\boldsymbol{k}]$は，要素の板厚を一定値$t$とし，応力-ひずみ関係行列を$[\boldsymbol{D_e}]$とすれば以下の通りとなる．

温度変化による節点力ベクトル

温度変化によるひずみ$\{\boldsymbol{\epsilon_0}\}$による節点力ベクトル$\{\boldsymbol{f_t}\}$は以下の通りとなる．

温度変化によるひずみの成分は以下のとおりである．

ここに，$\alpha$は線膨張係数，$\nu$はポアソン比，$T_p$は要素任意点の温度変化量，$\{\phi_i~~\phi_j~~\phi_k~~\phi_l\}^T$は節点温度変化量，$[N_1~~N_2~~N_3~~N_4]$は$[\boldsymbol{N}]$の要素である．また，温度変化量の符号は，温度上昇を正とする．

物体力による節点力ベクトル

要素内任意点に作用する物体力（ここでは慣性力とする）$\{\boldsymbol{F}\}=\{F_{Bx},F_{By}\}^T$は物体力に関わる節点力ベクトル$\{w\}$と内挿関数$[\boldsymbol{N}]$を用い以下の通り表せる．

ここに$\gamma$は要素の単位体積重量（質量$\times$重力加速度），$k_h$および$k_v$は水平方向・鉛直方向の慣性力の係数（重力加速度に対する比率）である．

以上より，物体力（慣性力）による等価節点外力ベクトル$\{\boldsymbol{f_b}\}$は，

として算定できる．ここでベクトル$\{\boldsymbol{w}\}$については$\gamma$をくくり出すととにより，

とおいている．

軸対称弾性応力解析

軸対称問題における応力-ひずみ関係とひずみ-変位関係

回転軸方向の応力・ひずみを添字$z$，半径方向の応力・ひずみを添字$r$，円周方向の応力・ひずみを添字$\theta$により表すとすれば，軸対称問題における応力-ひずみ関係式は，以下の通りとなる．

ここに，$\nu$はポアソン比，$\alpha$は熱膨張率，$T$は温度変化量である．

上式を変形することにより

これを，マトリックスおよびベクトルの形で表せば以下の通りとなる．

回転軸方向の変位および半径方向の変位を，それぞれ$w$，$u$とおけば，ひずみ-変位関係式は以下の通りとなる．

4節点4ガウス積分点パラメトリック要素としての定式化

ひずみ-変位関係行列$[\boldsymbol{B}]$の導入

軸対称要素内任意点の回転軸方向変位$w$および半径方向変位$u$を以下のように仮定する．ここに座標$(a,b)$は$[-1,1]$の範囲を有する正規化パラメータ座標，$w_{i,j,k,l}$および$u_{i,j,k,l}$は要素を構成する節点の変位を示す．

要素内任意点の変位を$\{\boldsymbol{u}\}$，節点変位を$\{\boldsymbol{u_{nd}}\}$とすれば，

要素内任意点のひずみ$\{\boldsymbol{\epsilon}\}$は，節点変位$\{\boldsymbol{u_{nd}}\}$を用いて，

要素内任意点の$r$座標値は，以下により評価する．

ここに，$r_i$，$r_j$，$r_k$，$r_l$は要素を構成する節点の$r$座標値である．

要素剛性行列

要素剛性行列は，本来要素の体積積分により求めますが，棒や平面問題では簡略化し1次元積分や2次元積分を計算することによりこれらを求めている．軸対称要素の場合を考えてみると，回転中心より半径$r$の位置にある微小要素の体積は，$r\times d\theta\times dr\times dz$となる．ここで積分値や荷重ベクトルなどを，円周方向は1ラジアンで評価することとし，Gauss積分を用いるための積分変数変換を行えば，

となる．よって，軸対称要素における要素剛性行列$[\boldsymbol{k}]$は，応力-ひずみ関係行列を$[\boldsymbol{D_e}]$とすれば以下の通りとなる．

２次元飽和-不飽和浸透流解析

２次元定常飽和浸透流問題の有限要素式

ダルシーの法則が成立する2次元直交異方性材料の直交座標系$x$-$y$における定常飽和浸透流問題の支配方程式は以下の通りである．なお奥行き方向の厚さは考慮しないことにする．

ここで，Galerkin法により，任意の$\delta \phi$に対し，支配方程式のの弱形式である次式を満足する未知数$\phi$を決定することを考える．

有限要素において，要素内の任意箇所の全水頭$\phi$が，内挿関数行列$[\boldsymbol{N}]$（ただし1行$\times$節点自由度数列）と節点全水頭ベクトル$\{\boldsymbol{h}\}$を用いて，

と表せれば，

となるが，$\{\boldsymbol{\delta h}\}$は任意の節点水頭ベクトルであるため，面積積分値が0となる必要がある．

次に部分積分公式

およびGreenの公式

を用いることにより，

また，$x$方向および$y$方向の流速をそれぞれ$v_x$，$v_y$とすれば，ダルシーの法則より

以上により，要素有限要素式は以下の通り表現できる．

なお$\{\boldsymbol{q}\}$は流速ベクトルの境界長さでの積分であり，流量を示すことは容易に理解できる．

4節点4ガウス積分点要素としての定式化

透水性マトリックス

面積積分の積分点を4点，線積分の積分点を2点とすれば，重みは1となるため，それぞれの行列・ベクトルは積分点での計算値の総和として次のように示される．

節点流量ベクトル$\{\boldsymbol{q}\}$については，等価節点ベクトルとして直接入力すればよい．

要素内平均流速

要素内任意点の全水頭$\phi$は内挿関数$[\boldsymbol{N}]$と節点全水頭ベクトル$\{\boldsymbol{h}\}$を用いて

であることから，ダルシーの法則より要素内任意点の流速$v_x$および$v_y$は

要素内平均流速は，要素内の4つのGauss積分点の平均とすればよいため，以下の通り求められる．

注意：以上は異方性地盤を対象とした展開であるが，プログラムでは，等方性地盤のみを扱える

不飽和浸透流解析への拡張

飽和浸透流解析

飽和定常浸透流解析は，以下の連立一次方程式を解く問題となる．

ここで，$\{q\}$は節点流量ベクトル，$\{h\}$は全水頭ベクトル，$[k]$は透水行列であり，以下の特性を有する．

• 透水行列$[k]$は定数，
• 節点流量ベクトル$\{q\}$の既知成分に相当する全水頭ベクトル$\{h\}$は未知数
• 節点流量ベクトル$\{q\}$の未知成分に相当する全水頭ベクトル$\{h\}$は既知数

飽和-不飽和浸透流解析

飽和-不飽和浸透流解析での連立方程式は，以下の特性を有する．

• 境界条件には，流量既知，全水頭既知，浸潤面が出現する境界の3種類を考慮する必要がある．
• 浸潤面が出現する境界以外は，流量か全水頭のいずれかが既知となる．
• 浸潤面が出現する境界上では流量・全水頭とも未知数となる．
• ただし浸潤面が出現する境界は，基本的に大気に接する境界であるため，雨の流入などを考慮しなければ流量は流出のみ(プログラム上は0以下)， 圧力水頭も0以下という制約がある．
• 飽和領域での透水行列の成分は定数であるが，不飽和領域での透水行列の成分はサクション圧(負の圧力水頭)の関数となっている．このため，収束計算により，未知数を定める必要がある．

収束計算の方法は，全節点に全水頭の初期値を与え，得られた解を次の収束計算の入力値とする，逐次代入法を用いている．この方法では全節点の水頭を仮定するため，仮定した水頭から全要素の透水行列を設定できるため，計算フローを単純化できるメリットがある．

全水頭の初期値は，全水頭既知境界以外は，飽和領域の最小全水頭を与える（要素領域の Z 座標の最小値）のが，解の精度・収束性を考慮したうえで有効のようである．

不飽和透水特性（van Genuchten モデル）

地盤の不飽和透水特性モデルは，飽和度-サクション圧関係，飽和度-不飽和透水係数比をテーブルで入力する方法もあるが，解析上，連続関数で取り扱うほうが，プログラムが単純化できるため，以下に示す van Genuchten モデルを採用している．

２次元非定常熱伝導解析

２次元非定常熱伝導問題の有限要素式

2次元等方性材料の直交座標系$x$-$y$における非定常熱伝導問題の支配方程式は以下の通りである．なお奥行き方向厚さは考慮しないことにする．

また，境界での熱流束$q$は，Fourierの法則より，境界上の外向き法線を$n$として，以下の通りとなる．

ここで，Galerkin法により，任意の$\delta T$に対し，支配方程式の弱形式である次式を満足する未知数$T$を決定することを考える．

有限要素において，要素内の任意箇所の温度$T$が，内挿関数行列$[\boldsymbol{N}]$（ただし1行$\times$節点自由度数列）と節点温度ベクトル$\{\boldsymbol{\phi}\}$を用いて，

と表せれば，

となるが，$\{\boldsymbol{\delta\phi}\}$は任意の節点温度ベクトルであるため，面積積分値が0となることが必要となる．

次に部分積分公式

およびGreenの公式

およびFourierの法則を用いることにより，$\kappa$に係る項は，

また

以上により，要素有限要素式は以下の通り表現できる．

境界条件に関する補足

境界（要素の辺）で熱流束$q_0$が与えられている場合は

とおけばよく，断熱境界では$q_0=0$とおけばよいことになる．

境界（要素の辺）で熱伝達がある場合は，熱伝達率を$\alpha_c$，境界温度を$T_c$として

となり，それぞれの成分を熱伝導マトリックスおよび熱流束ベクトルに足し込む必要がある．

境界条件を考慮した熱伝導マトリックス・熱流束ベクトルによる有限要素式

発熱体に関する補足説明

セメント・コンクリートをモデルとした発熱体を扱えます．セメント．コンクリートの断熱温度上昇量を，以下の式で仮定する．

上式の中の，$K$ (Tk) および $\alpha$ (Al) を材料特性として入力する．

上式を用い，発熱量 $Q$ および発熱率 $\dot{Q}$ は以下のように表現できる．

非定常有限要素式の解法

非定常有限要素式を解く上での時間離散化手法として，Crank-Nicolson差分式を用いる．

解析範囲全体の有限要素式を以下の通りとする．

ここで，時刻$t+\Delta t/2$での節点温度ベクトル$\{\boldsymbol{\Phi}\}$，節点温度ベクトル時間微分$\{\partial\boldsymbol{\Phi}/\partial t\}$，および節点熱流束ベクトル$\{\boldsymbol{F}\}$を以下の通りおきます．

これらを，上述の解析範囲全体の有限要素式に代入することにより，以下の差分式が得られ，時刻$t$における既知ベクトルを用いて$t+\Delta t$における節点温度を逐次求めることができる．

要素の材料特性は時間により変化しない設定をしているため，時間ステップ計算では，$numpy.linalg.inv(A)$ で一回だけ上式左辺の逆行列を求め，これを記憶して各時刻の温度履歴を計算している．

4節点4ガウス積分点要素としての定式化

面積積分の積分点を4点，線積分の積分点を2点とすれば，重みは1となるため，それぞれの行列・ベクトルは積分点での計算値の総和として次のように示される．なお，次式中$\ell$は境界条件を指定された要素辺の実長を表す．

要素熱伝導マトリックス

要素熱容量マトリックス

熱流束ベクトル

２次元弾性骨組構造振動解析

運動方程式とNermarkのβ法

単純化のため１質点系を考える．時刻$t+\Delta t$での質点系の釣り合いを示す運動方程式は以下の通り．ここにドットは時間微分を示す． \begin{equation} m\cdot\ddot{u}(t+\Delta t)+c\cdot\dot{u}(t+\Delta t)+k\cdot u(t+\Delta t)=f(t+\Delta t) \end{equation}

時刻$t+\Delta t$での変位$u$および速度$\dot{u}$をTaylor展開したものを以下に示す．

ここで，変位については右辺第4項において$1/3!=\beta$とおくことにより，速度については右辺第3項において$1/2!=\gamma$とおき，双方で高次項を無視することにより以下の式が得られる．

次に，$\dddot{u}$を以下のように線形化する．

上記関係を用いて，運動方程式を整理することにより，以下の方程式が得られる．

ここで，有限要素法への適用を考えると，変位を未知数として方程式を解きたいため，加速度と速度を既知の値および時刻$t+\Delta t$での変位で表すよう式を変形する．

上記を運動方程式に代入することにより，以下の関係が得られる．

ここで通常は，$\gamma=1/2$ が用いられ，$\beta=1/4$とすれば平均加速度法，$\beta=1/6$とすれば線形加速度法となる．

マトリックス表示した解くべき方程式

１質点系で求めた上式を，$\gamma=1/2$に固定し，マトリックス表示すると以下の通りとなる．

変位算定式

変位速度・加速度算定式

上記により変位$\{\boldsymbol{u}(t+\Delta t)\}$を求めた後，下記のより変位速度$\{\boldsymbol{\dot{u}}(t+\Delta t)\}$，加速度$\{\boldsymbol{\ddot{u}}(t+\Delta t)\}$を求める．

ここに，平均加速度法：$\beta=1/4$，線形加速度法：$\beta=1/6$である．

１質点系や多質点系においては，加速度を未知数として解く式をよく見かけるが，有限要素法への適用においては，非線形を取り入れる場合に有利であるなどの理由から，変位を未知数とする場合も多いようである．

また直感的に釣り合い式の表現は，$m\cdot a=f$か，$k\cdot u=f$ですが，筆者の直感的感覚は，有限要素法では$k\cdot u=f$であり，プログラムは変位を未知数として解くよう作成することにした．

剛性マトリックス・質量マトリックス

座標変換マトリックス$[\boldsymbol{T}]$は，全体座標系変位を要素座標系変位に変換するものである． 全体座標系での剛性マトリックス$[\boldsymbol{K}]$と要素座標系での剛性マトリックス$[\boldsymbol{k}]$には次の関係がある．

全体座標系での質量マトリックス$[\boldsymbol{M}]$と要素座標系での質量マトリックス$[\boldsymbol{m}]$についても同様に次の関係がある．

減衰マトリックス

実際の計算で，構造物の減衰を考慮する場合，減衰マトリックス$[\boldsymbol{C}]$を具体的に与える必要がある． ここでは，減衰マトリックスの与え方について考えてみる．

単純化のため地盤加速度を受ける１質点系の運動方程式を考える．

ここで，よく使われる減衰の形として，Rayleigh減衰を考える．これは，減衰係数は質量と剛性に比例するという考え方であり，以下の式表現となる．

上式を考慮して減衰定数$h$を書き直してみると ($\omega$は固有円振動数，$f$は振動数)

よって，もし現地測定などで，$i$次モードと$j$次モードの減衰定数$h_i$，$h_j$および固有振動数$f_i$，$f_j$が判れば，上式を連立させることにより，比例定数$\zeta_m$および$\zeta_k$を求めることができる．

減衰係数をマトリックス表示すれば，上記で定めたスカラー量$\zeta_m$および$\zeta_k$を用いて，

なお，実務上は比例定数の求め方として，以下のような手法が用いられる場合がある．

• 振動モードに応じた減衰定数がよくわからないので，減衰定数の値は1つとして，一般に使われている値（例えば$h_1=h_2=0.05$）などにセットする．
• 固有振動数については，固有値解析などで求めることができるため，1次および2次の固有振動数$f_1$，$f_2$を求める．
• 2つの固有振動数$f_1$，$f_2$と2つの減衰定数$h_1=h_2=0.05$が定まれば，上式により，比例定数$\zeta_m$および$\zeta_k$を求めることができる．

固有値解析

構造系の固有円振動数 $\omega$ をもとめるには，一般化固有値問題として，以下の固有方程式を解けば良い．

固有円振動数 $\omega_n$ が求まれば，以下の関係により固有周期 $f_n$ を求めることができる．

２次元弾性骨組構造幾何学的非線形解析

増分形剛性方程式

よく知られているように，幾何学的非線形項を含む増分形剛性方程式は，以下の通りとなる．

内力計算

内力計算

内力増分 $\{\Delta f^*\}$ は，線形剛性マトリックス $[K_L]$ 及び剛体回転を除去した部材座標系での変位増分 $\{\Delta u^*\}$ から計算される． \begin{equation*} \{\Delta f^*\}=[K_L]\{\Delta u^*\} \end{equation*} \begin{align*} &\{\Delta f^*\} =\begin{Bmatrix} \Delta N^*_i & \Delta S^*_i & \Delta M^*_i & \Delta N^*_j & \Delta S^*_j & \Delta M^*_j \end{Bmatrix}^T \\ &\{\Delta u^*\}=\begin{Bmatrix} \Delta u^*_i & \Delta v^*_i & \Delta \theta^*_i & \Delta u^*_j & \Delta v^*_j & \Delta \theta^*_j \end{Bmatrix}^T \end{align*} \begin{align*} &\Delta u^*_i=0 &~& \Delta v^*_i=0 &~& \Delta \theta^*_i=(\tan\theta^*_i)_k-(\tan\theta^*_i)_{k-1} \\ &\Delta u^*_j=\Delta\ell & & \Delta v^*_j=0 & & \Delta \theta^*_j=(\tan\theta^*_j)_k-(\tan\theta^*_j)_{k-1} \end{align*}

ここに， $\Delta\ell$ は収束計算中の前回と今回の要素長の変化を表す． 回転量についても，前回回転量（添字 $k-1$）と今回回転量（添字 $k$）が考慮される．

剛体回転の除去

剛体回転の除去方法を以下に示す．

正接加法定理を用いて，

Elimination of rigid rotation

上式を用いて，現時点での $(\tan\theta_i^*)_k$, $(\tan\theta_j^*)_k$ および前回の $(\tan\theta_i^*)_{k-1}$, $(\tan\theta_j^*)_{k-1}$ を計算できる． これらから，以下のように，剛体回転を除去した回転量増分が計算できる．

これらのステップでは，初期状態（変形前）の節点座標から計算した座標変換行列を用いる．

弧長法

図に示すような，単純化されたスカラー量の荷重-変位曲線を考える． この図を参照して，以下の関係が導かれる．

Concept of Arc-Length method

$\Delta\lambda$ の初期値

$\Delta U_R=0$ を仮定すれば，$\Delta\lambda$ は以下のように表せる．

上式は，プラス・マイナスの２つの値を取り得る．また弧長法ではこの$\Delta\lambda_0$の符号が非常に重要になる．

$\Delta\lambda_0$ の符号は，以下のように定める．

• 前回釣り合い点から今回釣り合い点までの変位増分ベクトル ${\Delta U_{-1}}$ を定義する
• $\{\Delta U_0\}$ の計算（$\{\Delta U_0\}=[K_T]^{-1}\{\Delta F\}$）
• 内積 $\{\Delta U_{-1}\}^T\{\Delta U_0\}=|\Delta U_{-1}|\cdot |\Delta U_0|\cdot \cos\theta$ の計算．ここで， $\theta$ は２つの変位ベクトルのなす角を示す．
• 内積が，$\{\Delta U_{-1}\}^T\{\Delta U_0\} \geqq 0$ なら，角度 $\theta$ は90度以下である．この場合，$\Delta\lambda_0$ の符号は正（プラス）とする．
• 内積が，$\{\Delta U_{-1}\}^T\{\Delta U_0\} \lt 0$ なら，角度 $\theta$ は90度より大きい．この場合，$\Delta\lambda_0$ の符号は負（マイナス）とする．

スケーリングパラメータ $\phi$ は，以下の式で計算する．このプログラムでは， $\alpha=1.0$ としている． \begin{equation*} \phi=\sqrt{\frac{\alpha}{\{\Delta F\}^T\{\Delta F\}}} \end{equation*}

弧長 $\Delta s$ は，$\Delta\lambda=1.0$ を仮定して以下の式で求める．

収束計算のための補正係数 $\Delta\lambda$

前出の図を参照して， $\Delta s$ を $\Delta L$ に置き換えることにより，以下の式が得られる．

上式の $(\Delta L)^2$ を最小化する条件より，$\Delta\lambda$ は以下のように計算できる．

解析のフローチャート

Flowchart of Calculation

実際のプログラムでは， $\Delta s$ と $\Delta\lambda$ は以下のような扱いとなっている．

• 初期載荷重（nnn=1）において， $\Delta s$ を計算し，また $\Delta\lambda_0=1.0$ にセットする．
• もし計算された $\Delta\lambda$ が $\Delta\lambda_0$ と比較して非常に小さい場合， $\Delta s$ は増加させる必要がある．なぜならば $\Delta s$ は $\Delta\lambda=1.0$ を仮定して計算されているからである．
• もし計算された $\Delta\lambda$ が初期の $\Delta\lambda_0 (=1.0)$ より大きければ，$\Delta s$ の増加はストップされる．なぜならば過大な $\Delta s$ は，予期せぬ挙動を引き起こす原因になるからである．

# Initial parameter setting for Arc-length method
if nnn==1:
    ds0=np.sqrt(np.sum(dis0*dis0)+spara*spara*np.sum(df*df))
    ds=ds0
    dlam0=1.0
dlam=np.sqrt(ds*ds/(np.sum(dis0*dis0)+spara*spara*np.sum(df*df)))
if np.abs(dlam) < 0.1*np.abs(dlam0): ds=ds*1.2
if np.abs(dlam0) < np.abs(dlam): ds=ds
dlam=np.sqrt(ds*ds/(np.sum(dis0*dis0)+spara*spara*np.sum(df*df)))
if np.sum(dis_ref*dis0)<0.0: dlam=-dlam

２次元弾性骨組構造座屈解析

固有方程式

座屈荷重および変形モードを求めるための固有方程式は，以下のとおりである．

$[K_L]$および$[K_G]$は，幾何学的非線形解析に用いるものと同じ形であるが，$[K_G]$に含まれる軸力$P$について，座屈解析では圧縮を正とするところに注意が必要である．

下記固有方程式を解くことにより係数$\lambda$が求まるが，$[K_G]$に含まれる軸力$P$を確定する必要があるため，単位の荷重を入力して，固有値解析の前に，線形解析を１回実施しておく必要がある．

１次元非定常熱伝導解析

非定常有限要素式の解法

２次元非定常熱伝導解析と同様，Crank-Nicolson 差分式に有限要素式を当てはめた結果は，以下のとおりである．

それぞれのマトリックス・ベクトルを，要素単位で示すと以下の通りとなる．

発熱体に関する説明

セメント・コンクリートをモデルとした発熱体を扱える．セメント．コンクリートの断熱温度上昇量を，以下の式で仮定している．

上式の中の，$K$ (Tk) および $\alpha$ (Al) を材料特性として入力する．

上式を用い，発熱量 $Q$ および発熱率 $\dot{Q}$ は以下のように表現できる．