最新版を正確に読む場合には pdf ファイル をどうぞ。これは web 検索のための簡易旧版です。

• 5.3.1 要素剛性方程式
  • 5.3.1.1 柱の場合
  • 5.3.1.2 梁の場合
  • 5.3.1.3 等価節点外力の意味と値
  
  
• 5.3.2 直接剛性法と内力分布
  • 5.3.2.1 直接剛性法
  • 5.3.2.2 2要素による例題
  • 5.3.2.3 曲げモーメント等の内力分布
  
  
• 5.3.3 平面トラス
  • 5.3.3.1 全体座標系と要素座標系
  • 5.3.3.2 不静定トラスの例題
  
  
• 5.3.4 平面骨組
  • 5.3.4.1 全体座標系での要素剛性方程式
  • 5.3.4.2 2部材の不静定骨組の例題
  
  
• 5.3.5 全体剛性方程式の解法
• 5.3.6 その他の要素
  • 5.3.6.1 弾性床上の梁の要素
  • 5.3.6.2 梁の動的問題に対する要素

5.3 構造部材の有限要素

5.3.1 要素剛性方程式

5.3.1.1 柱の場合

図 5.8: 柱の1有限要素

前節の区分的多項式による近似を一般化し，図-5.8の ような長さ$\ell$の1有限要素に外力が作用して 図のような変位が生じてつり合ったとして近似方程式を求めてみよう。 この一般的な場合の仮想仕事の式 は

$$\int_0^\ell EAu' \delta u'\dint x -\int_0^\ell p \delta u \dint x -F_1 \delta u_1-F_2 \delta u_2 = 0$$ (5.13)

となる。ただし，ここの$\ell$は1有限要素の長さである。 つまり前節のような区分的多項式による近似の1要素に注目し， 相隣接する要素同士はあとで適切に連続させることにしているだけで， 本質的には前節での定式化との違いは無い。 まず，この1有限要素の変位を多項式で近似することから始めよう。 変断面でも容易に定式化できるが， 簡単のために1有限要素は単一材料でできた等断面部分だとする。

式(5.13)で，変形できる物体にとって 最も重要な項は第1項の内力仮想仕事項だから，それが 無くならないためには，$u(x)$の近似関数には1階の微係数が存在する必要がある。 したがって，多項式を用いる限りは`$x$'の1次以上の 高次の多項式が必要となる。また境界条件によっては 柱が剛体的に移動する場合もあるため，平行移動を表す`$1$'（定数）も必要である。 さらに幾何学的境界条件からは， 隣接する要素との間で変位$u(x)$が両端で連続する必要があるから，1有限要素には 最低でも自由度が二つ必要になる。 このように考えると，最も次数の低い多項式の選択としては， 二つの自由度$c_1$, $c_2$を用いて

$$u(x)\sim c_1+c_2 x$$ (5.14)

と置くことが可能なようだ。ここでも やはり，物理的に意味の無い係数$c_n$を使う代わりに，前節で示したように 両端の変位$u_1$, $u_2$を未知係数とするような表現にしたい。 そうすることによって，幾何学的境界条件を満足させられる近似関数と 対応する重み関数にすることが容易にできそうだ。つまり上式から

$$u_1=u(0)=c_1, \qquad u_2=u(\ell)=c_1+c_2 \ell$$

という関係を得るので，結局$c_1=u_1$, $c_2=\left(u_2-u_1\right)/\ell$と 置き換えればよく，式(5.14)は次のように書き改められる。

$$u(x)\sim u_1\phi_1(x)+u_2\phi_2(x), \qquad \phi_1(x)\equiv \dfrac{\ell-x}{\ell}, \quad \phi_2(x)\equiv \dfrac{x}{\ell}$$ (5.15)

このような関数$\phi_n(x)$は変位を近似するために仮定した 関数なので変位関数 ，あるいは一般的には試行関数 とか内挿関数 と^5.3呼ばれる。

以上を仮想仕事の式(5.13)に代入して整理すると

$$\int_0^\ell EA\left(u_1\phi_1'+u_2\phi_2'\right) \left(\de... ...u_2\phi_2\right) \dint x -\delta u_1 F_1-\delta u_2 F_2 = 0$$

となる。すなわち各仮想変位$\delta u_1$と$\delta u_2$毎にまとめると

$$\delta u_1\left\{ l_{11}u_1+l_{12}u_2-p_1-F_1 \right\}+ \delta u_2\left\{ l_{21}u_1+l_{22}u_2-p_2-F_2 \right\}=0 \eqno{(*)}$$

を満足していればいいことになる。ここに

$$p_i \equiv \int_0^\ell \phi_i(x) p(x)\dint x, \quad l_{ij}\equiv \int_0^\ell EA \phi_i'(x) \phi_j'(x) \dint x$$ (5.16)

で定義したものはそれぞれ，等価節点外力 成分および要素剛性行列 要素である。後者は具体的にその積分を実行すると

$$l_{11}=l_{22}=\dfrac{EA}{\ell}, \quad l_{12}=l_{21}=-\dfrac{EA}{\ell}$$ (5.17)

という値を持つ。したがって 式($*$)が任意の仮想変位$\delta u_1$と$\delta u_2$に対して 成立するためには，それぞれの係数がいずれも零であればよく

$$\left\{\begin{array}{r}F_1 F_2\end{array}\right\}+ \left\{... ...ay}\right)\left\{ \begin{array}{r}u_1 u_2\end{array}\right\}$$ (5.18)

を得る。この式は，一つの要素の両端節点に作用している力と節点変位の 関係を代数方程式で表したつり合い式であり，要素剛性方程式と 呼ばれる。特にこの式は引張り圧縮のみで抵抗する柱の1要素の剛性方程式である ため，柱の要素剛性方程式 と呼ばれる。式(5.17) (5.18)から 柱の1有限要素は，バネ定数 $(\slfrac{EA}{\ell})$の線形バネとして 取り扱えばいいこともわかる。 だから，章-2でトラス部材のことを検討したときに， 簡単のために線形バネに置き換えたのである。

5.3.1.2 梁の場合

図 5.9: 梁の1有限要素

梁の場合も図-5.9のような 一般的な状態を対象として仮想仕事式 を求めると

$$\int_0^\ell EIw'' \delta w''\dint x -\int_0^\ell q \delta w \dint x$$ (5.19)

$$\mbox{}\qquad -S_1 \delta w_1-C_1 \delta\theta_1 -S_2 \delta w_2-C_2 \delta\theta_2 = 0$$

となる。ただし$\ell$は1有限要素長である。$w(x)$を近似する関数は

1. 1要素の剛体的な沈下も表現できるためには，`1'が必要となる。
2. また1要素の剛体的な（微小）回転も表現できるためには，`$x$'が必要。
3. 式(5.19)で最も重要な内力仮想仕事項が 無くならないためには，2階の微係数が存在する必要があり，`$x^2$'以上が 必要となる。
4. 要素間あるいは端部でたわみとたわみ角が連続あるいは規定 できるためには，一つの要素は四つの自由度を持つ必要があり， たわみを近似するときの未知係数としては，両節点のたわみとたわみ角の 合計四つを用いればいい。

という条件から，四つの自由度を持つ3次の多項式が最低次の多項式であり

$$w(x)\sim a_1+a_2x+a_3x^2+a_4x^3$$

と近似すればよさそうだ。 さらに，この物理的に意味の無い未知係数$a_i$ ($i=1\sim 4$)を 要素両節点でのたわみとたわみ角に関係付けると，最終的に

$$w(x)\sim w_1\psi_1(x)+\theta_1\psi_2(x)+w_2\psi_3(x)+\theta_2\psi_4(x)$$ (5.20)

という表現になる。 ここに梁の変位関数は

$$\left\{\begin{array}{ll} \psi_1(x)\equiv 1-\dfrac{3x^2}{\ell^2} +... ...3}, & \psi_4(x)\equiv \dfrac{x^2}{\ell} -\dfrac{x^3}{\ell^2} \end{array}\right.$$ (5.21)

と定義される。

式(5.20)を式(5.19)に代入して整理すると

$$\begin{eqnarray*} &&\delta w_1 \left\{ ( k_{11} w_1+k_{12} \theta_1+k_{13}\... ...k_{42} \theta_1+k_{43} w_2+k_{44} \theta_2)-q_4-C_2\right\}=0 \end{eqnarray*}$$

となる。ここに$k_{ij}$と$q_i$は

(5.22)

で定義される剛性行列要素と等価節点外力 成分である。これが任意の$\delta w_n$, $\delta\theta_n$ ($n=1, 2$)に 対して成立する条件から，梁の要素剛性方程式 が

$$\left\{\begin{array}{c}S_1 C_1 S_2 C_2\end{array}\right\... ...t) \left\{\begin{array}{c}w_1 \theta_1 w_2 \theta_2\end{array}\right\}$$ (5.23)

と表現できる。

図 5.10: 梁のたわみに対して仮定した区分的多項式

式(5.21)を式(5.22a)に代入すると， 剛性行列$\fat{k}_b$ が具体的に

$$\fat{k}_b\equiv\mat{k_{ij}}= EI\left(\begin{array}{cccc} \... ...icolumn{3}{l}{\mbox{Symm.}}&\dfrac{4}{\ell} \end{array}\right)$$ (5.24)

と求められる。 なお簡単のためにこの章では，括弧無しの太字で行列を表している。 また，変位関数を図-5.4と同様の区分的多項式として 表したのが図-5.10である。

5.3.1.3 等価節点外力の意味と値

図 5.11: 等価節点外力の意味

一つの例として，両端固定の棒に図-5.11のような 等分布外力のみが作用している場合に，端部の反力を 剛性方程式(5.18) (5.23)で計算してみる。 この場合は剛性方程式右辺の両端の変位成分はすべて零になるので， 柱と梁の支点反力はいずれも

$$\left\{\begin{array}{r}F_1 F_2\end{array}\right\}=- \left\... ... q_3 q_4\end{array}\right\} \renewedcommand{arraystretch}{1}$$

となる。すなわち等価節点外力の符号を換えたものは それぞれ，分布外力を受ける両端固定の棒の不静定反力であることがわかる。 特に等分布外力である場合の値を計算してみると

$$p_1=p_2=\dfrac{p_0\ell}{2}, \quad q_1=q_3=\dfrac{q_0\ell}{2}, \quad q_2=-q_4=-\dfrac{q_0\ell^2}{12}$$ (5.25)

という値になる。$p_1$, $p_2$, $q_1$, $q_3$は分布外力の長さ方向の総和の 半分ずつになっているので直感的な理解も容易である。 また$q_2$, $q_4$は節-4.1.4 (2)で 算定されていた不静定モーメント反力式(4.33)に一致して おり，1要素剛性方程式から求められる両端の不静定反力も厳密解に 一致していることがわかる。 ただし，このように反力が厳密解に一致するのは，節点変位が 厳密解に一致するのと同様，梁や柱特有の特殊な場合のみに 見られることであり， 一般には近似解しか得ることができないことには十分留意しておく必要がある。

1. 演習問題5-3の6番で 求めた境界値問題の仮想仕事式に基づいて$u(x)$に 仮定すべき試行関数の候補を挙げ，それを用いて剛性行列と 等価節点外力を求めよ。

5.3.2 直接剛性法と内力分布

5.3.2.1 直接剛性法

前節では，対象とする構造系を複数の有限要素に分割し，あとで適切に 連続させることを条件とした上で個々の要素で仮想仕事式が成立する ものとして，1有限要素のつり合い式に対する要素剛性方程式を導いた。 ここでは，その分割した要素同士を適切に連続させる方法と， それを用いて元の境界値問題を解く方法とを説明し，いくつかの例を示す。 実際には，剛性行列や等価節点外力も電子計算機の中で数値的に 計算されれることも多く，以下の手法も数値的に処理されることが多いが， ここではその基礎について梁を対象とした具体例を用いて説明する。

図 5.12: 節点での連続条件と直接剛性法

いま図-5.12のように二つの有限要素の共有 節点2に集中外力が作用している場合を考えて，この要素を 連続させて全体剛性方程式を求める方法を説明しよう。図のように この節点2のすぐ左右をそれぞれ`$2-$', `$2+$'という 節点番号で区別すると，それぞれの要素毎に式(5.23)の 要素剛性方程式が成立するので

$$\left\{\begin{array}{c}S_1 C_1 S_{2-} C_{2-}\end{array... ... \theta_1 w_{2-} \theta_{2-}\end{array}\right\} \eqno{(a)}$$

$$\left\{\begin{array}{c}S_{2+} C_{2+} S_3 C_3\end{array... ... \theta_{2+} w_3 \theta_3\end{array}\right\} \eqno{(b)}$$

を満足している。ここに上付きの(l) (r)はそれぞれ，節点2のすぐ 左右の要素の値を区別するために用いた記号である。 載荷された部分では図の下に示したように，力とモーメントのつり合いから

$$S_2=S_{2-}+S_{2+}, \qquad C_2=C_{2-}+C_{2+} \eqno{(c)}$$

を満足するように，左右の要素で外力を分担して抵抗している。一方， この部分のたわみとたわみ角がもし連続していないと，二つの 要素がここで分離したり折れ曲がったりしてしまう。したがって 変位の連続条件は

$$w_{2-}=w_{2+}=w_2, \qquad \theta_{2-}=\theta_{2+}=\theta_2$$

でなければならない。 変位の連続条件を考慮しながら式(a)の第3, 4式と 式(b)の第1, 2式とを 上の連続条件式(c)に代入して整理すると

$$\left\{\begin{array}{c}S_2 C_2\end{array}\right\}= \left(\... ... \theta_3\end{array}\right\} \renewedcommand{arraystretch}{1}$$

となる。これと式(a)の上2行と，式(b)の下2行の 四つの式とを併せて行列表示すると

$$\left\{\begin{array}{c}S_1 C_1 S_2 C_2 S_3 C_3\end... ... \theta_1 w_2 \theta_2 w_3 \theta_3\end{array}\right\}$$

が二つの要素を繋いだ全体剛性方程式になっている。

もう少し見通しがよくなるように次の記号を定義する。

$$\fat{f} \equiv \left\{\begin{array}{c} S C \end{array}\rig... ...persc{b})\supersc{t} & \fat{k}_b\supersc{c} \end{array}\right)$$ (5.26)

ここに， $\fat{k}_b\supersc{x}$ (X=A, B, C)は それぞれ剛性行列$\fat{k}_b$の2$\times$2の小行列である。 これを用いて上の最後の結論を記述すると，二つの要素で全体系を 解こうとする場合の全体剛性方程式は

$$\left\{\begin{array}{c} \fat{f}_1 \fat{f}_2 \fat{f}_3... ...y}{c} \fat{w}_1 \fat{w}_2 \fat{w}_3 \end{array}\right\}$$ (5.27)

となる。すなわち下線部のように，相隣接する要素の共有節点での 各量を加え合わせれば，最終的な全体剛性方程式が求められることが わかる。このように行列を機械的に重ねていく方法を直接剛性法 と呼んでいる。

5.3.2.2 2要素による例題

図 5.13: 2要素の場合の直接剛性法の例

具体的な使い方を図-5.13の例を用いて説明しよう。 この場合は境界条件から，まず

$$w_1=0, \quad \theta_1=0, \quad w_3=0, \quad \theta_3=0$$

であり，外力条件は

$$S_2=Q, \quad C_2=0$$

となる。これを式(5.27)に代入すると，その 第1, 2, 5, 6式は左辺右辺共に未知量を含んでいるのに 対し，第3, 4行目のみは左辺がすべて既知であり，未知量は 右辺にしか存在しないことがわかる。よって後者の2式を 取り出すと

$$\left\{\begin{array}{c}Q 0\end{array}\right\}= \left(\begi... ...ight) \left\{\begin{array}{c}w_2 \theta_2\end{array}\right\}$$

を得る。ここでは左右の要素の剛性と長さが等しいので， それを区別する上付き記号は省略した。 具体的に値を代入して逆行列を計算すると

$$\left\{\begin{array}{c}w_2 \theta_2\end{array}\right\}= \f... ...{\begin{array}{c}\slfrac{Q\ell^3}{24EI} 0\end{array}\right\}$$

と中央のたわみとたわみ角を求めることができる。 この解を元の全体剛性方程式の第1, 2, 5, 6式に 代入することによって，支点反力が $S_1=S_3=\slfrac{-Q}{2}$, $C_1=C_3=\slfrac{Q\ell}{4}$と求められ，たわみ・たわみ角も支点反力・不静定 モーメントも，厳密解（図-4.29付近）の$\ell$を$2\ell$に 置き換えた答に一致する。

5.3.2.3 曲げモーメント等の内力分布

図 5.14: 連続梁を三つの有限要素で解く例

もう一つの簡単な例として図-5.14にある 等分布外力を受ける2径間連続梁を解いてみる。 特徴を出すために，長さの違う左右のスパンの曲げ剛性が 異なる値を持つようにした。図では梁の上面を揃えてあるが， 正確には2部材の図心の位置も揃っている^5.4ものとする。 式(5.25)で等分布外力の 場合の等価節点外力成分は算定できているので， 式(5.26)を用いてそれぞれの要素の 要素剛性方程式を簡易表示すると

$$\left\{\begin{array}{c} \fat{f}_1 \fat{f}_{2-}\end{array}\... ...ay}\right\},                                       \eqno{(d)}$$

$$\left\{\begin{array}{c} \fat{f}_{2+} \fat{f}_{3-}\end{arr... ...ft\{\begin{array}{c} \fat{w}_3 \fat{w}_4\end{array}\right\}$$

でなければならない。 また境界条件および載荷条件は

$$w_1=w_2=w_4=0, \quad \theta_1=\theta_4=0, \quad S_3=0, \quad C_2=C_3=0$$

で与えられる。

上の三つの要素剛性方程式を用いて直接剛性法を使い， それに上で示した境界条件および載荷条件を代入すると， 全体剛性方程式が

$$\left\{\begin{array}{@{}c@{}} S_1 C_1 S_2 0 0 0... ...a_2 w_3 \theta_3\\ 0 0 \end{array}\right\} \eqno{(e)}$$

となる。この第4, 5, 6式だけが右辺のみに 未知量を含む式になっており，それを取り出すと

$$\begin{eqnarray*} \left\{\begin{array}{c}0 q_0\ell 0\end{array}\right\}&=& \... ...ft\{\begin{array}{c}\theta_2 w_3 \theta_3\end{array}\right\} \end{eqnarray*}$$

となっている。したがって逆行列を必要な成分だけ計算すると最終的に変位が

$$\left\{\begin{array}{c}\theta_2 w_3 \theta_3\end{array}\... ... \slfrac{37}{2016} \slfrac{1}{(224\ell)}\end{array}\right\}$$

と求められる。これを残りの全体剛性方程式($e$)に代入すれば， 支点反力と不静定モーメントが

$$S_1=-\dfrac{2q_0\ell}{7}, \quad C_1=\dfrac{q_0\ell^2}{84}, \... ..._4=-\dfrac{121q_0\ell}{112}, \quad C_4=-\dfrac{65q_0\ell}{168}$$

と求められる。

求められた節点変位の値を式(5.20)で仮定した 変位関数の係数に代入することによって，各要素の任意点での 変位を関数として得ることができるから，その2階の微係数を算定すれば 各要素毎の曲げモーメント分布も求めることができる。 しかし，変位関数が3次の多項式で近似されていたため，曲げモーメント 分布は線形にしかならず，この例の場合も計算すると 左側のスパンで

$$M=-2EIw''=\dfrac{q_0\ell}{14} \left(2\ell-6x\right)$$

となり，右側のスパンも二つの要素のそれぞれの左端を$x$の原点にすれば

$$M=\dfrac{q_0\ell}{112}\left(-16\ell+47x\right), \quad M=\dfrac{q_0\ell}{112}\left(31\ell-65x\right)$$

にしかならない。図-5.15の直線分布が このようにして計算されたモーメント分布であり， 曲線が厳密解を示している。図-5.6の柱の場合には 軸力が最小2乗法 的近似解になっていたが，梁の場合には有限要素法の 曲げモーメントの近似解，つまりたわみの2階の微係数が， 厳密解の最小2乗法的な近似解になっていることがこの図から明らかである。

図 5.15: 曲げモーメント図の求め方

ところで，有限要素法が近似解法であるにもかかわらず， 柱と梁の場合は，節点で得ることができた変位成分は厳密解に一致していたし， 等分布外力の作用した両端固定梁の不静定反力も厳密解に一致していた。 つまり上で求めたように，式($e$)の全体剛性方程式からは 支点反力と不静定モーメントの正解が求められていた。 それなら，同じように曲げモーメント図も厳密なものが計算できないだろうか。 そこで今度は要素剛性方程式($d$)に注目すると，この要素剛性方程式の 両端の外力は要素の内力そのもののか，あるいはその符号を替えたものに 一致していることが，その定式化や図-5.12等からわかる。 そこで，求められた変位の解を上述のように全体剛性方程式や変位関数の微係数に代入する代わりに，要素剛性 方程式($d$)に代入してみよう。 そして内力の正の向きと外力の向きとを 考慮しながら，各要素毎に端部の曲げモーメントを要素剛性方程式から 計算すると

$$-C_1=-\dfrac{q_0\ell^2}{84}, \quad C_{2-}=-C_{2+}=-\dfrac{19q_0\ell^2}{84}, \quad C_4=-\dfrac{65q_0\ell^2}{168}$$

と求められる。実は$C_1$と$C_4$は上で求めた厳密な不静定曲げモーメントに 一致しており，また$C_{2-}$と$-C_{2+}$は中間節点箇所の左右の 内力曲げモーメント（この例では集中外力モーメントが作用していないから 同じ値）になっている。 そして，静定系の曲げモーメント分布が（有限要素法を使うまでもなく） 既知であることから，一旦節点の内力曲げモーメントが得られさえすれば， 簡単に図-5.15の曲線のような 厳密な曲げーモーメント分布図を描くことができるはずだ。 つまり梁の問題では，正しい曲げモーメント図を計算するためにも 有限要素法は有用な手法の一つであると言えよう。

このように，有限要素法は近似解法の一つではあるものの， 対象が柱か梁に限定されていれば，節点での変位成分と要素剛性方程式から 求められる内力の値とは厳密解に等しいことを示すことができた。 現在のように有限要素法が一般的な境界値問題の汎用的近似解法と して認識される以前に，節-5.4.1に示すように， 元の境界値問題を厳密に解いて剛性方程式と同じ形式の基礎式を求め， それを用いた解法をマトリックス構造解析法 と呼んで構造力学・構造工学分野では用いていた。 あるいは古典的には節-5.5に示す たわみ角法や三連モーメントの定理として用いていた。 この基礎式と有限要素法の剛性方程式とは完全に一致するか， 単なる項の並べ替えを行ったものになっていることから， 有限要素法で節点での厳密解が求められるのは当然なのである。 有限要素法の発達はこのマトリックス構造解析法を発端としている。 ただ一般論として，有限要素法があくまでも近似解法であることは 心に留めておくべきである。

図 5.16: 数値計算を用いた例

次の演習問題のようなまっすぐの連続梁を解く有限要素プログラムは 各種存在すると思われるが，著者らが作成したプログラムも電子的に 入手可能になっている。図-5.16のような 等分割要素に対する処理をパーソナルコンピュータ 上で会話型で行えるようになっている。 スパン途中にヒンジやバネ支持も挿入できる。 入手方法についてはまえがきを参照のこと。 できれば各自，プログラミングの練習として少なくとも一直線上の連続梁くらいは 解けるものを作って欲しい。 このような連続梁であれば，直接剛性法をプログラミングすることは そんなに難しいものにはならないからだ。 また，有限要素法を越えて，上で説明したように曲げーモーメントと せん断力の分布図を描くデータも計算するように工夫してみて欲しい。

図 5.17: 不静定系の有限要素解析の練習

1. 先端に集中外力が作用した片持ち梁を一つの有限要素で 解き，厳密解と比較せよ。
2. 図-5.17の(a)の変断面両端固定の 梁を二つの有限要素で解き，節点2のたわみ・たわみ角を 計算した上で，要素剛性方程式を利用して曲げモーメント図を求めよ。
3. 図-5.17の(b)の連続梁を三つの 有限要素で解き，節点2のたわみ・たわみ角，節点3のたわみ角を 計算し，全体剛性方程式から支点反力を求めよ。結果を要素剛性方程式に 代入して正しい曲げモーメント図を求めよ。
4. 同じく図-5.17の(b)の連続梁で， 節点3の支点が$\Delta$だけ沈下する場合を三つの 有限要素で解き，未知変位と支点反力および曲げモーメント図を求めよ。

5.3.3 平面トラス

5.3.3.1 全体座標系と要素座標系

図 5.18: 平面不静定トラスの例
図 5.19: トラス部材の有限要素

柱が解けるのなら図-5.18の不静定平面トラスも 有限要素法で解けるだろう。 トラスの1本1本の部材は両端がヒンジで接続されており，部材中間には 外力が作用していないため，曲げモーメントによってではなく 軸の圧縮引張りだけで抵抗している。したがって トラスの各部材を柱要素として取り扱うことができるはずだ。

さて， 図に示したようにトラスの各部材はお互いに異なる方向を向いているため， 単純な柱と違い部材両端にはせん断外力成分も作用させることができる。 しかし上でも述べたように，このせん断外力に対して曲げで 抵抗するとは考えない上に， 部材に分布する外力も考えない。 したがって，図-5.19のように それぞれの成分を定義したとしても，トラス部材の要素剛性方程式は 柱としての抵抗しか無いので

$$\left\{\begin{array}{c} F_1 S_1 F_2 S_2 \end{array}\ri... ...ft\{\begin{array}{c} u_1 w_1 u_2 w_2 \end{array}\right\}$$ (5.28)

と書かざるを得ない。つまり，剛性行列の第2, 4行列の要素は零 であり，両端自由の場合にはせん断外力には抵抗できない。

図 5.20: 平面トラス有限要素の要素座標系と全体座標系

一方それぞれの部材は同じ方向を向いたものではなく， それぞれが任意の方向を向いている。 つまり，このトラス部材の場合の要素剛性方程式(5.28)は それぞれの部材軸方向を要素座標とする系で定義されたものである。 ということは，ある一つの節点での力と変位成分が， 隣り合う部材間では異なる方向の成分として定義されてしまうことになる。 このままでは，複数の有限要素同士の連続条件を考慮するのに，座標系の 違いを考慮して計算する等煩雑な手続きが必要となる。 そこで，すべての要素に共通する全体座標系での成分を各節点の 未知数とするような要素剛性 方程式を誘導しておけば，前節で用いた直接剛性法が，お互いに向きの 異なる要素間に対してもそのまま適用でき便利である。

すなわち，図-5.20の$\xi$-$\zeta$軸で定義される 要素座標系においては式(5.28)が成立している。 一方，$x$-$z$座標系は空間に固定した座標系で， すべての要素に共通する全体座標系とする。ということは， 要素剛性方程式(5.28)の両辺の各量を$x$-$z$全体座標系で定義された 量に変換すれば，すべての部材の要素剛性方程式を，一つの共通した 全体座標系における量で表現できることになる。 柱の場合はたまたま，剛性方程式両辺の列行列が ちょうど2次元ユークリッド空間のベクトル成分と 同じになっているから，この二つの座標の成分間の関係は

$$\left\{\begin{array}{l} \overline{u}_n=u_n\cos\alpha+w_n\si... ...overline{S}_n=-F_n\sin\alpha+S_n\cos\alpha \end{array} \right.$$

となる$(n=1,2)$。したがって，座標変換行列を

$$\fat{T}_2\equiv\left(\begin{array}{cc} \cos\alpha&-\sin\alpha \sin\alpha&\multicolumn{1}{r}{\cos\alpha} \end{array}\right)$$ (5.29)

と定義しておくと，それぞれのベクトル成分が

$$\left\{\begin{array}{c} u_1 w_1 u_2 w_2 \end{array}\ri... ...t\{ \begin{array}{c} F_1 S_1 F_2 S_2 \end{array}\right\}$$

となる。この両式を式(5.28)の両辺に代入して 整理すると，最終的な全体座標系における要素剛性方程式 が

$$\left\{\begin{array}{c} \overline{F}_1 \overline{S}_1\\ ... ...e{w}_1\\ \overline{u}_2 \overline{w}_2 \end{array}\right\}$$ (5.30)

と求められる。 あるいはこの右辺の行列乗算を実行し，座標変換された剛性行列 としてその成分を陽に表すと

$$\dfrac{EA}{\ell}\left(\begin{array}{cccc} \cos^2\alpha&-\si... ...ulticolumn{3}{l}{\mbox{Symm.}}&\sin^2\alpha \end{array}\right)$$ (5.31)

となっている。

5.3.3.2 不静定トラスの例題

例として図-5.18の平面トラスを解いてみよう。 部材1-2は長さが $\slfrac{\sqrt{3}\ell}{2}$であるが， 要素座標系が全体座標系と一致しているので$\alpha=0$でいい。 一方部材2-3は要素座標系を節点2から3の方向に とった方が便利なので，全体座標系とのなす 角度は $\alpha=-150^\circ$あるいは$210^\circ$になっている。 したがって，それぞれの要素剛性方程式は式(5.31)を 用いて

$$\left\{\begin{array}{c}\overline{F}_1 \overline{S}_1\\ \... ...{w}_{2}\\ \overline{u}_3 \overline{w}_3\end{array}\right\}$$

と表現できる。この二式ともに同じ全体座標系で定義されているから 共有節点2に対して直接剛性法を用いることができ，最終的な 全体剛性方程式は

$$\left\{\begin{array}{c}\overline{F}_1 \overline{S}_1\\ \... ...e{w}_2\\ \overline{u}_3 \overline{w}_3 \end{array}\right\}$$

となる。下線部が共有節点2について重ね合わせた部分である。 境界条件は，節点1, 3ですべての変位成分が固定されていること， 節点2に二つの外力が作用していることから

$$\overline{u}_1=\overline{w}_1=\overline{u}_3=\overline{w}_3=0, \quad \overline{F}_2=P, \quad \overline{S}_2=Q$$

と与えられている。全体剛性方程式のうち， 左辺が与えられている第3, 4行目が先に解けるはずなので

$$\left\{\begin{array}{c}P Q\end{array}\right\}= \dfrac{EA}{... ...{array}{c} \overline{u}_2 \overline{w}_2 \end{array}\right\}$$

となり，これを節点2の変位成分について解けばいい。その結果

$$\left\{\begin{array}{c}\overline{u}_2 \overline{w}_2\end{a... ...{\sqrt{3}P}{4}+\slfrac{(8\sqrt{3}+9)Q}{12} \end{array}\right\}$$

を得る。これを残りの全体剛性方程式の第1, 2, 5, 6行目に代入 すれば，支点反力がそれぞれ

$$\left\{\begin{array}{c}\overline{F}_1 \overline{S}_1\end{a... ...t\}= \left\{\begin{array}{c} \sqrt{3}Q -Q\end{array}\right\}$$

と算定できる。

設計のことを考えると， 最終的にはそれぞれの部材の軸力を計算する必要があるが，これも 連続梁の例と同様， 要素剛性方程式の$\xi$方向の両端の節点軸力$F_{2-}$, $F_3$， あるいはその符号を替えたもの$-F_1$, $-F_{2+}$で算定できる。 あるいは上で得た全体座標系での変位成分を要素剛性方程式に代入して 求められる要素両端外力に，座標変換行列を逆にかけて 要素座標系に戻すことによっても節点軸力を求めることができる。 この例では，要素1-2は $N=-\overline{F}_1=P+\sqrt{3}Q$となる。 要素2-3では座標変換が必要になり

$$N=-\overline{F}_3\cos 30^\circ+\overline{S}_3\sin 30^\circ=-2Q$$

となる。図-5.21に結果を示した。

図 5.21: 2本トラスの変位と軸力

図 5.22: 不静定平面トラス-- 1パネル長10m，高さ8m

図-5.22に示したのは，パーソナルコンピュータ 上で作った有限要素プログラム による連続不静定トラスの解析例で ある。a〜oの15本の部材が節点1〜9で繋がった 簡単なトラスで，各部材の長さと外力は図示した通りだが，Young率は すべて200GN/m$\mbox{}^2$と共通にした。ただし 断面積は表-5.1の3種類を 用いている。図-5.22の 上図には全体剛性方程式を解いて得た変位を誇張して 示し，下の図には軸力を棒グラフで示した。 また，求められた軸力と支点反力の数値も表-5.2に 列挙してある。負の軸力は圧縮で，上弦材が圧縮を受けている。 反力は$z$方向の外力として示したので，すべて上向きの反力である。 このプログラムもまえがきに書いた方法で入手可能である。 例として図-5.22用のデータファイル等を添付してある。

表 5.1: 例の部材断面積 (m$^2$)

部材	断面積
a, d, h, l, o	0.01
b, f, j, n	0.008
c, e, g, i, k, m	0.004

表 5.2: 例の部材軸力と反力 (kN)

部材	軸力	部材	軸力	部材	軸力	部材	軸力	節点	反力
a	$-$26	e	$33$	i	$-62$	m	$56$	1	$-22$
b	14	f	$10$	j	$26$	n	$29$	5	$-81$
c	26	g	$-33$	k	$62$	o	$-56$	9	$-47$
d	$-$28	h	$7$	l	$-59$	--	--	--	--

1. 図-5.23の左の平面トラスを三つの有限要素を 用いて解け。要素座標の向きを表すときの角度の正方向に注意せよ。 また，三つの部材が集まる点での直接剛性法について熟慮せよ。
2. 同じ図の右側に示した5本トラスを，プログラムを用いて 数値的に解いてみよ。ただし$E=200$GN/m$\mbox{}^2$, $A=0.005$m$\mbox{}^2$と する。

図 2.9: 平面トラス

5.3.4 平面骨組

5.3.4.1 全体座標系での要素剛性方程式

骨組というのは，梁と柱を縦横に並べ，それぞれを剛結して 組み立てた構造系である。 さらに，ある平面内に組み立てられた骨組で，その同じ面内にのみ外力が 作用して，その面内でのみ変形する骨組の問題を平面骨組問題と呼ぶ。 したがって各部材は，軸力の作用と曲げの作用を同時に受けると 考えなければならず，上で誘導した柱および梁の要素剛性方程式を同時に 考慮する必要がある。 そこで図-5.24のような1有限要素を考えることにすると， 式(5.18)の柱の剛性方程式と式(5.23)の梁の 剛性方程式とを一本化して

$$\left\{\begin{array}{c}F_1 S_1 C_1 F_2 S_2 C_2\end... ...1 w_1 \theta_1\\ u_2 w_2 \theta_2\end{array}\right\}$$ (5.32)

とすればいいことがわかる。剛性行列の各要素および等価節点外力は 各成分共に既に定義されている。 等価節点外力については，与えられている量であって 演算上は節点の集中外力と同様に取り扱えばいいから，以下，簡単の ために省略する。

図 5.24: 平面骨組の有限要素
図 5.25: 平面骨組有限要素の要素座標系と全体座標系

骨組の場合も一般には前節の平面トラスと同様，各部材はそれぞれ別々の 方向を向いて配置されている。トラスとの違いは，部材間が たいていの場合は剛結されている点である。ヒンジが入るような 場合も簡単に定式化できるが，詳細については 文献[118]等を参照のこと。 式(5.32)は図-5.25の$\xi$-$\zeta$座標系 で成立する要素剛性方程式である。要素座標系と全体座標系の 間の関係の中でトラスに無かったのは曲げモーメントとたわみ角で あるが，ここで対象としているのは平面問題なのでこの二つに ついては座標変換の必要が無い。したがって，座標変換行列を

$$\fat{T}\equiv\left(\begin{array}{ccc} \cos\alpha & -\sin\al... ...icolumn{1}{r}{\cos\alpha} & 0\\ 0 & 0 & 1 \end{array}\right)$$ (5.33)

と拡張して定義すると，それぞれの 節点外力および変位の座標変換則は

$$\left\{\begin{array}{c} \overline{u}_i \overline{w}_i\\ ... ...}} \left\{\begin{array}{c} F_i S_i C_i \end{array}\right\}$$ (5.34)

となる$(i=1,2)$。 したがって，平面トラスの場合と同様の座標変換の演算を 行うことによって，最終的な全体座標系における要素剛性方程式 が

$$\left\{\begin{array}{c} \overline{\fat{F}}_1 \overline{\... ...verline{\fat{u}}_1 \overline{\fat{u}}_2 \end{array}\right\}$$ (5.35)

と表現できる。 ここに，全体座標系での外力ベクトル成分と変位成分を簡単のために

$$\overline{\fat{F}}_i\equiv \lfloor \overline{F}_i\;\; \over... ...}_i\;\; \overline{\theta}_i \rfloor\supersc{t} \quad (i=1,2)$$ (5.36)

と置いた。$\fat{k}_f$は式(5.32)右辺にある骨組の剛性行列 で

$$\fat{k}_f\equiv \renewedcommand{arraystretch}{0.7} \left(\be... ...persc{b})\supersc{t} & \fat{k}_f\supersc{c} \end{array}\right)$$ (5.37)

と定義した。 最後の式の $\fat{k}_f\supersc{x}$ (X=A, B, C)は 実線で区分けした3$\times$3の小行列である。 この小行列を用いると式(5.35)は

$$\left\{\begin{array}{c} \overline{\fat{F}}_1 \overline{\... ...verline{\fat{u}}_1 \overline{\fat{u}}_2 \end{array}\right\}$$ (5.38)

とも表現できる。

5.3.4.2 2部材の不静定骨組の例題

図 5.26: 平面骨組の例

図-5.26の例を解いてみよう。 この場合は要素1-2は$-30^\circ$傾いた部材であるが，要素2-3は 全体座標系と要素座標系が一致している。したがって， 直接剛性法によって全体剛性方程式を求めると

$$\left\{\begin{array}{@{}c@{}} \overline{\fat{F}}_1 \over... ...overline{\fat{u}}_2 \overline{\fat{u}}_3 \end{array}\right\}$$

でいい。ここの座標変換行列は

$$\fat{T}_1\equiv \fat{T}(\alpha=-30^\circ)=\left(\begin{array... ...} & \slfrac{\sqrt{3}}{2} & 0 \\ 0 & 0 & 1 \end{array}\right)$$

という値を持つ。 境界条件と載荷条件は

$$\overline{u}_1=\overline{w}_1= \overline{u}_3=\overline{w}_... ...verline{F}_2=0, \quad \overline{S}_2=Q, \quad \overline{C}_2=0$$

と与えられているので，上の全体剛性方程式のうち左辺が既知の第4, 5, 6行目が 先に解くべき式になり

$$\left\{\begin{array}{c} 0 Q 0 \end{array}\right\}= \Mat{... ...2 \overline{w}_2\\ \overline{\theta}_2 \end{array}\right\}$$

を解けばいいことになる。具体的に剛性行列の要素に値を代入すると

$$\left\{\begin{array}{c} \overline{u}_2 \overline{w}_2\\ ... ...ht)^{-1} \left\{\begin{array}{c} 0 Q 0 \end{array}\right\}$$

となるから，この逆行列を計算すれば解を得る。 その結果を用いて全体あるいは要素剛性方程式を前節同様 適切に用いると，支点反力や曲げモーメント図・せん断力図・ 軸力図が算定できる。

5.3.5 全体剛性方程式の解法

ここまでの例で見たように，剛性方程式を用いて解く場合にはまず 力学境界条件の与えられている節点での式が解け，すべての節点変位が 先に求められる。よく変位法 と呼ばれる理由はここにある。 ここでは，計算機の中で剛性行列に境界条件・外力条件を 代入して変位を計算する手順について概説しておく。

全体座標系における集中外力ベクトルと分布外力ベクトルが それぞれ $\overline{\fat{f}}^c$, $\overline{\fat{f}}^d$で 表され，直接剛性法によって全体剛性方程式を

$$\overline{f}_i\equiv\overline{f}_i^c+\overline{f}_i^d =\sum... ...ine{k}_{ij} \overline{u}_j, \qquad i=1,2,\cdots N \eqno{(a)}$$

のように得たとする。これに外力条件を代入するのは簡単で， 左辺の $\overline{\fat{f}}$の対応する成分を与えればいい。 次に，幾何学的境界条件が例えば節点$n$で

$$\overline{u}_n=\Delta_n \eqno{(b)}$$

と与えられた場合，これを上の全体剛性方程式(a)に 代入して整理すると

$$\overline{f}_i-\overline{k}_{in}\Delta_n= \sum_{j=1}^{n-1}\... ...+1}^N \overline{k}_{ij} \overline{u}_j, \quad i=1,2,\cdots N$$

となる。この第$n$行目のみを式(b)と置き換え， その他の行はそのままにして行列表示すると

$$\left\{\begin{array}{@{}c@{}} \overline{f}_1 \vdots\\ \... ...\!+\!1} \vdots \vdots \overline{u}_N \end{array}\right\}$$ (5.39)

を全体剛性方程式と考えればいいことになる。 第$n$行目は式(b)そのものである。 このようにして， 与えられた幾何学的境界条件をすべて全体剛性方程式に組み込んでしまえば， 左辺は既知量のみとなり， 正しい組み合わせで幾何学的境界条件が与えられてさえいれば 右辺の剛性行列も正則になって逆行列が存在し，すべての変位成分を 求めることができる。このような手順で全体剛性方程式を解き， 内力を要素剛性方程式から求める場合の代表的な フローチャートを図-5.27に示した。 また図-5.28(a)は， 数値的に解析した骨組の例である。 この解析プログラム もまえがきに記した方法で入手可能である。 与えられる境界条件の種類に制限があるが，支点沈下等も取り扱うことが できるようになっている。

図 5.27: 有限要素解析のフローチャート例

図 5.28: 平面骨組の解析例題

1. 図-5.28(b)の正三角形の平面骨組を，三つの有限要素を 用いて解け。直接剛性法において，剛性行列をどのように 重ねていけばいいか注意深く考えよ。
2. 図-5.28(a)の平面骨組を， プログラムを用いて数値的に解いてみよ。 ただしYoung率は$E=200$GN/m$\mbox{}^2$とし， 部材2-6, 4-7を除いて$A=0.001$m$\mbox{}^2$, $I=0.0001$m$\mbox{}^4$とし，部材2-6, 4-7の面積と 断面2次モーメントは共にこの2倍とする。

5.3.6 その他の要素

5.3.6.1 弾性床上の梁の要素

節-4.7.1で取り上げた弾性床上の梁についても， 同様の手順で要素剛性方程式を誘導できる。仮想仕事式 が

$$\int_0^\ell EIw'' \delta w''\dint x + \int_0^\ell k_w w \d... ..._1 \delta\theta_1 -S_2 \delta w_2-C_2 \delta\theta_2 = 0$$ (5.40)

であることは，梁の場合と同様の演算を式(4.86)に 対して行えば容易に示すことができる。 用いる変位関数は，その条件が初等梁理論と同じであることから 前節で用いた変位関数式(5.20)が使え， 最終的に剛性方程式 が

$$\renewedcommand{arraystretch}{0.7} \left\{\begin{array}{c}S_... ... \theta_2\end{array}\right\} \renewedcommand{arraystretch}{1}$$ (5.41)

と求められる。ここに$\fat{k}_b$は式(5.24)で 定義されたものと同じであり，上式(5.40)左辺第2項に 関する剛性行列が

$$\fat{k}_w \equiv \left( \left[\int_0^\ell k_w \psi_i(x) \... ...mn{3}{l}{\mbox{Symm.}}&\slfrac{\ell^3}{105} \end{array}\right)$$ (5.42)

で定義される。これが，$k_w$が一様な場合の弾性床抵抗を代表する剛性行列 である。 弾性床上の梁の問題の変位の厳密解に指数関数が存在することは， 式(4.87)で示した通りである。 したがって，区分的多項式で変位を近似した 剛性方程式から求められる有限要素法による近似解は，上の柱や梁の場合とは異なり， 節点位置であっても厳密解には一致しないことに注意すべきである。 解法は普通の梁と同じなので割愛する。

図 5.29: 弾性床上梁の解の収束

図-5.29には， 弾性床上の単純梁の中央に集中荷重を載せた場合の， 中央のたわみ$w$と曲げモーメントの有限要素解と 要素数$N$の関係を示した。 曲げモーメントは，変位関数の2階の微係数で求めた 場合 $\left(-EIw''\right)$と，図-5.15で示したのと同様に 要素剛性方程式の反力から求めた場合 $\left(M\right)$との二つを示した。 縦軸は有限要素解と厳密解の相対誤差 $\left(\mbox{Err}\sub{r}\right)$であり， 両軸ともに対数表示である。 要素サイズを小さくする（要素数を増やす）に従って，近似誤差 が急速に零に漸近している。 図中の直角三角形の勾配が$2$であるから，変位関数の微係数で求めた 曲げモーメント $\left(-EIw''\right)$が 要素分割数の2乗の速度で零に収束していることがわかる。 また，図-5.7, 5.15の場合と同様に， たわみの高階微係数で求めた断面力の誤差の方がたわみの誤差よりも大きく， 収束の速度も遅い。 ただし，要素剛性方程式から求めた断面力 $\left(M\right)$の誤差と 収束速度は，変位そのものの誤差の傾向とほぼ同様になっているのは興味深い。

5.3.6.2 梁の動的問題に対する要素

梁の運動方程式は式(4.88)なので， この初期値境界値問題に対しても静的問題と 同様に弱形式を求めてみると，形式的に式(5.55)のように なることは容易にわかる（正確な誘導は節-10.4.7を 参照のこと）。 ここで解くべき問題は時間と場所の関数を求めることであるから， この「動的な仮想仕事式」に対する動的な剛性方程式を 評価するための変位の近似では，式(5.20)を少し変更して

$$w(x,t)\sim w_1(t) \psi_1(x)+\theta_1(t) \psi_2(x) +w_2(t) \psi_3(x)+\theta_2(t) \psi_4(x)$$ (5.43)

と置いてみればよさそうだ。 ちょうど変数分離で初期値境界値問題を解くのに似ている。 以下，粘性項は無いものとして誘導する。 式(5.43)を式(5.55)に代入して演算すると， 静的な剛性方程式と異なる項は慣性項のみである。 これは

$$\begin{eqnarray*} m \int_0^\ell \ddot{w} \delta w \dint x &=& m \bigl[ \d... ...ta_1 + \mu_{43} \ddot w_2 + \mu_{44} \ddot \theta_2 \} \bigr] \end{eqnarray*}$$

となる。ここに

$$\mu_{ij} \equiv \int_0^\ell \psi_i(x) \psi_j(x) \dint x$$ (5.44)

と定義した。 したがって要素の運動方程式 は

$$\fat{m}_b \ddot{\fat{w}}_b(t) + \fat{k}_b \fat{w}_b(t) = \fat{f}_b(t)$$ (5.45)

となる。 ここに$\fat{m}_b$は $\left(m \mu_{ij}\right)$を$ij$要素とする質量 行列は，具体的には式(5.42)の$k_w$を$m$で置き換えて

$$\fat{m}_b \equiv \left( \left[\int_0^\ell m \psi_i(x) \ps... ...mn{3}{l}{\mbox{Symm.}}&\slfrac{\ell^3}{105} \end{array}\right)$$ (5.46)

という成分を持ち，$\fat{k}_b$は式(5.24)で定義した 梁の曲げに関する剛性行列である。また

$$\fat{w}_b(t) \equiv \lfloor w_1(t) \; \theta_1(t) \; w_2(t)... ...+ \lfloor q_1(t)\; q_2(t)\; q_3(t)\; q_4(t) \rfloor\supersc{t}$$

と定義した。

式(5.45)は1要素の動的な要素剛性（運動）方程式であり， 境界条件を代入する前なら 形式的に4自由度系の運動方程式と同じ形をしている。 もちろん4質点系のそれ（対角行列）とは異なり， 質量行列は式(5.42)のように非対角項にも非零な値を 持つ4$\times$4の行列になっている。 このような質量行列を整合質量行列 と呼ぶ。梁の動的問題をこの要素剛性方程式を用いて解く場合には， 質量行列についても剛性行列同様の直接剛性法を適用し， 共有節点に関する要素同士を重ね合わせてやればいい。

自由振動の問題を例として解いておこう。この場合には 外力を零にしておいて周期解を探せばいいだけなので

$$\fat{w}_b(t)=\fat{a} \exp(\omega t)$$

と置けばいい。ここに$\fat{a}$は未定の振幅列行列である。 これを式(5.45)の要素剛性方程式に代入すると， 形式的には

$$\left( \fat{k}_b - \omega^2 \fat{m}_b \right) \fat{a} = \fat{0}$$

を満足していればいいことになる。したがって， 有意な$\fat{a}$が存在するための条件から

$$\det \vert \fat{k}_b - \omega^2 \fat{m}_b \vert = 0$$ (5.47)

でなければならず，この式が固有振動数$\omega$を 決定する特性方程式になる。 上で「形式的」と書いたのは，簡単のために境界条件等の 処理についての記述を割愛しているからである。 複数の有限要素で解く場合には，直接剛性法でそれぞれの 行列を重ねたあとに境界条件を代入し，そのあとで上の演算を 行う必要がある。処理法の詳細については座屈についての同様の 計算(p.)あるいは節-10.4.7を 参照のこと。

表 5.3: 固有振動数の近似度の改善

要素数	$\zeta_1$	$\zeta_2$	$\zeta_3$	$\zeta_4$	$\zeta_5$
1	1.10992	1.27157	--	--	--
2	1.00395	1.10992	1.23994	1.27157	--
4	1.00026	1.00395	1.01827	1.10992	1.12909
8	1.00002	1.00026	1.00129	1.00395	1.00927
16	1.00000	1.00002	1.00008	1.00026	1.00063

このようにすると，式(5.47)のように， 行列の固有値問題 の固有値として系の固有振動数を計算できる。 両端単純支持梁の固有振動数を，用いる要素数を増やしながら 求めたのが表-5.3である。 ただし，表で用いた

$$\zeta_k \equiv \dfrac{\mbox{有限要素解}}{\mbox{厳密解}}= \dfrac{\omega_k}{(k\pi/\ell)^2 \sqrt{EI/m}}$$

は，第$k$次振動数の厳密解と近似解の比である。 要素数が増えるにつれて 解の精度が上がっていくことが明らかである。 工学的に最も重要な最低次の固有振動数は，1要素（2自由度）であっても わずか11%程度の誤差しか含んでおらず，2要素（4自由度）で は1%未満になっていることは興味深い。 振動数の近似解が厳密解より大きくなるのは，近似することによって モデルの変形が真の変形より拘束されたものになってその剛性が大きく なってしまうからである。

より複雑な系の場合には，上の整合質量行列の代わりにそれぞれの節点に 何らかの集中質量を配置した集中質量行列 を算定して用いることが，特に自由振動解析において行われることがある。 これについても別途参考文献を参照のこと。 また，自由振動問題ではなく強制振動問題の時刻歴応答解を求める場合には， 式(5.45)を時々刻々解いていかなければならず， 時間微分の項をある時間刻みで差分表示のような漸化式に置き換えて 数値計算することになる。数値解析上の安定性や精度等の問題点や， 具体的な解法等については節-10.4.7や 参考文献[8]等を参照のこと。

最新版を正確に読むためには pdf ファイル をどうぞ。これは web 検索のための簡易旧版です。