を入れて,,, について解けば ( ) ( ) 得る. よって となるが ( / ( ) ( ) と無次元化している ), これを { N ( ) } { d} と表現して内部変位 と節点変位 { d} とを結びつける { } { ( ) ( ) } (.) (.) 節点での F と M は図. の

日本機械学会宇宙工学部門

第３章 有限要素法（その１）

３.１ 梁要素

有限要素法においては外力も境界条件も節点で考える．もちろん分布荷重は考慮でき るが，要素上の分布荷重は適当に節点への等価は集中荷重として置き換える．こう考える と梁の曲げの方程式(2.8)において分布荷重無し(p(x)=0)の

0

dx

w

d

4 4

=

EI

(3.1) である．この一般解は先に解いたように 3 3 2 2 1 0

c

x

c

x

c

x

c

w

=

+

+

+

(3.2) となる．梁の有限要素の長さを L とすると，その両端 x=0,L にて境界条件を導入して解い てきた． いまここで考えるのは梁の中の一部をきりとった要素である．未知数が４つ(

c

₀ ∼

c

₃) であるので両端で２つづつの境界条件を入れればよい． 入れ方として

M

or

0

dx

dw =

F

or

w

=

0

の２つの量の組合わせが考えられるが，右と左の節点で採用する物理量が異なっては統一 性がとれない．従ってひとつの方法として w と dw/dx の変位を採用する方法があり，こ れを変位法という．さて，この境界条件 1

w

w

=

,

₁

dx

dw

b

=

at

x

=

0

,

w

w

=

₂

=

b

₂

dx

dw

at

x

=

L

日本機械学会宇宙工学部門

w

1

b

1

w

b

∂ w_{∂ x} を入れて

c

₀

,

c

₁

,

c

₂

,

c

₃について解けば 1 0

w

c

=

1 1

c

=

b

(

2 1

)

(

1 2

)

2 2

2

L

1

w

w

L

3

c

=

-

-

b

+

b

(

1 2

)

2

(

1 2

)

3 3

L

1

w

w

L

2

c

=

-

+

b

+

b

得る．よって

ï

ï

þ

ï

ï

ý

ü

ï

ï

î

ï

ï

í

ì

b

b

´

2 2 1 1

w

w

（3.3） となるが(

x

=

x

/

L

と無次元化している)，これを

( )

{

}

{ }

d

=

N

x

w

（3.4） と表現して内部変位 w と節点変位

{ }

d

とを結びつける

{ }

N

を変位関数と呼ぶ． 節点での F と M は図 3.1 の符号と図 2.3 とを見比べて

( )

3 3 3 1

6

c

dx

w

d

0

F

F

=

-

=

EI

=

EI

( )

3 2

F

L

6

c

F

=

=

-

EI

( )

2 2 2 1

2

c

dx

w

d

0

M

M

=

=

-

EI

=

-( )

L

2

c

6

c

L

M

M

2

=

-

=

2

+

3

(

)

(

)

{

2 3 2 3 2 3 2 3

}

2

3

2

2

3

1

-

x

+

x

x

-

x

+

x

x

-

x

-

x

+

x

=

L

L

w

日本機械学会宇宙工学部門

であるので，節点での値を添字 1,2 で表すと

ï

ï

þ

ï

ï

ý

ü

ï

ï

î

ï

ï

í

ì

b

b

ú

ú

ú

ú

û

ù

ê

ê

ê

ê

ë

é

-EI

=

ï

ï

þ

ï

ï

ý

ü

ï

ï

î

ï

ï

í

ì

2 2 1 1 2 2 2 2 3 2 2 1 1

4

6

2

6

6

12

6

12

2

6

4

6

6

12

6

12

w

w

L

L

L

L

L

L

L

L

L

L

L

L

L

M

F

M

F

(3.5) となる．これを

{ }

F

=

[ ]

K

{ }

x

(3.6) と表示して剛性方程式と呼び，[K]を剛性行列という． 対 称 行 列 に な っ て い る こ と に 留 意 さ れ た い ． こ れ は 偶 然 で な く ， 物 理 的 に は Maxwell-Bettiの相反定理によるものであり，3.4 節のエネルギー原理による導出によ り明らかになるであろう．一定の力 F に対して K が大きければ変形 x は小さくなるので K が大きいほど剛な構造ということができる．分布荷重がなければこれが厳密解である．従 って，分布荷重のない一様な梁の場合，いくら長くても要素分割は１つでよい．分布荷重 がある場合，一つの梁をいくつかの要素に分割して分布荷重を等価な集中荷重に置き換え る．これについては次の例題で実例にそって考える．

3.1.1 例１ 片持梁

L

P

F =-P

2

b

1

=0

=0

左の固定点において

w

₁

=

b

₁

=

0

，右の自由端で

F

₂

=

-

P

,

M

₂

=

0

であるの で剛性方程式は

日本機械学会宇宙工学部門

ï

ï

þ

ï

ï

ý

ü

ï

ï

î

ï

ï

í

ì

b

ú

ú

ú

ú

û

ù

ê

ê

ê

ê

ë

é

-EI

=

ï

ï

þ

ï

ï

ý

ü

ï

ï

î

ï

ï

í

ì

-2 2 2 2 2 2 3 1 1

0

0

4

6

2

6

6

12

6

12

2

6

4

6

6

12

6

12

0

w

L

L

L

L

L

L

L

L

L

L

L

L

L

P

M

F

(3.7) である．この式を計算すると上半分と下半分について剛性行列が

þ

ý

ü

î

í

ì

b

ú

û

ù

ê

ë

é

-EI

=

þ

ý

ü

î

í

ì

2 2 3 1 1

2

6

6

12

w

L

L

L

L

M

F

(3.7)

þ

ý

ü

î

í

ì

b

ú

û

ù

ê

ë

é

-EI

=

þ

ý

ü

î

í

ì-2 2 2 3

4

6

6

12

0

w

L

L

L

L

P

(3.8) となる．上半分の式(3.8)は変位

w

₂

,

b

₂ が得られたあと反力を計算する式で，下半分の 式(3.9)は外力が与えられたときに変位

w

₂

,

b

₂を計算する式である．この後者(3.9)に ついて未知数

w

₂

,

b

₂ について解くと

þ

ý

ü

î

í

ì-ú

û

ù

ê

ë

é

-÷

ø

ö

ç

è

æ EI

=

þ

ý

ü

î

í

ì

b

-0

P

L

4

L

6

L

6

12

L

w

1 2 1 3 2 2

þ

ý

ü

î

í

ì

-=

þ

ý

ü

î

í

ì-ú

û

ù

ê

ë

é

EI

=

EI EI 2 3 3 3 2 1 2 1 2 1 3 1 3

0

PL PL L L L

P

L

となる．また節点１での反力

F

₁

,

M

₁ は得られた節点変位

w

₂

,

b

₂ を剛性方程式の右上半 分の式(3.8)に適用して

þ

ý

ü

î

í

ì

=

þ

ý

ü

î

í

ì

-ú

û

ù

ê

ë

é

-EI

=

þ

ý

ü

î

í

ì

EI EI

PL

P

L

2

L

6

L

6

12

L

M

F

2 PL 3 PL 2 3 1 1 2 3 となる．これらの結果は前章の例１の厳密解(2.15)と一致する．

3.１.2 例２:両端単純支持

この例では中央に集中荷重がある．有限要素法では力と変位は節点においてのみ与えら れるから中央点にも節点を設け，２要素３節点でモデル化することになる．節点１と２， 節 点 ２ と ３ と で 作 ら れ る 剛 性 行 列 を 重 ね 合 せ , 力 の 条 件

M

1

=

M

2

=

M

3

=

0

,

P

F

₂

=

-

と変位の条件

w

₁

= w

₃

=

0

を代入すると

日本機械学会宇宙工学部門

w

₁

₌₀

_{w =0}

3

P

図 3.3 ：両端単純支持梁

L

2

=



ú

ú

ú

ú

ú

ú

ú

ú

û

ù

ê

ê

ê

ê

ê

ê

ê

ê

ë

é

-+

+

-+

-+

-EI

=

ï

ï

ï

ï

þ

ïï

ï

ï

ý

ü

ï

ï

ï

ï

î

ïï

ï

ï

í

ì

2 2 2 2 2 2 2 2 3 3 2 1

L

4

L

6

L

2

L

6

L

6

12

L

6

12

L

2

L

6

L

4

L

4

L

6

L

6

L

2

L

6

L

6

12

L

6

L

6

12

12

L

6

12

L

2

L

6

L

4

L

6

L

6

12

L

6

12

L

0

F

0

F

0

F

ï

ï

ï

ï

þ

ïï

ï

ï

ý

ü

ï

ï

ï

ï

î

ïï

ï

ï

í

ì

b

b

b

´

3 2 2 1

0

w

0

を得る．ここで節点２において両要素が重なっていることに注意されたい．例１の式(3.9) と同じく未知節点変位に関して第 2,3,4,6 番目の行と列を抜出して剛性方程式を再構成 すると

ï

ï

þ

ï

ï

ý

ü

ï

ï

î

ï

ï

í

ì

b

b

b

ú

ú

ú

ú

ú

û

ù

ê

ê

ê

ê

ê

ë

é

-EI

=

ï

ï

þ

ï

ï

ý

ü

ï

ï

î

ï

ï

í

ì

-3 2 2 1 2 2 2 2 2 2 2 3

w

L

4

L

2

L

6

0

L

2

L

8

0

L

2

L

6

0

24

L

6

0

L

2

L

6

L

4

L

0

0

P

0

日本機械学会宇宙工学部門

この式を変位について解けばよいが，文字計算は面倒なので数値例として



=2 (したがっ て L=



/2=1) を採用して計算すると

ï

ï

þ

ï

ï

ý

ü

ï

ï

î

ï

ï

í

ì

ú

ú

ú

ú

û

ù

ê

ê

ê

ê

ë

é

-EI

-=

ï

ï

þ

ï

ï

ý

ü

ï

ï

î

ï

ï

í

ì

b

b

b

-0

0

1

0

4

2

6

0

2

8

0

2

6

0

24

6

0

2

6

4

P

w

1 3 2 2 1

ï

ï

þ

ï

ï

ý

ü

ï

ï

î

ï

ï

í

ì

-EI

-=

4 1 6 1 4 1

0

p

を得る．前章の同じ例題の結果(2.18)で



=2 とおけば

EI

-=

EI

-=

=

6

P

48

P

w

w

3 max 2



となって厳密解に一致する． また

b

₂=0 となって対称性も満たしている．ここでは表示 上，剛性方程式を解くのに逆行列を使っているが，実際には連立１次方程式を解けばよい．

3.1.3 例３:分布荷重のかかる片持ち梁

分布荷重がかかっているときにはモデル化するわけにはいかない．なぜなれば有限要素 を導出するときに採用した変位関数は分布荷重がないとして導かれたものだからである． よってこの例題の場合，分布荷重を節点上の集中荷重に置き換えてやらなければならない． この例題の場合，４つの要素に分割して，分布荷重の合計 p



をその 1/4 づつの要素の両

p

L

4



日本機械学会宇宙工学部門

端に振り分けて図のように等価な力をかけることにする．



=4 として数値計算を行う．先ほどの２つの例と同様に未知変位に関する行と列を残して

ï

ï

ï

ï

ï

þ

ï

ï

ï

ï

ï

ý

ü

ï

ï

ï

ï

ï

î

ï

ï

ï

ï

ï

í

ì

b

b

b

b

ú

ú

ú

ú

ú

ú

ú

ú

ú

ú

ú

û

ù

ê

ê

ê

ê

ê

ê

ê

ê

ê

ê

ê

ë

é

-EI

-=

ï

ï

ï

ï

ï

þ

ï

ï

ï

ï

ï

ý

ü

ï

ï

ï

ï

ï

î

ï

ï

ï

ï

ï

í

ì

5 5 4 4 3 3 2 2

w

w

w

w

4

6

12

2

6

8

.

sym

6

12

0

24

2

6

8

6

12

0

24

2

6

8

6

12

0

24

0

1

0

2

0

2

0

2

8

p

が得られる．この剛性行列をみればわかるように対角の近傍を除いてゼロである疎行列と なっている．実用的な有限要素法プログラムはこの行列の性質を使って記憶容量と演算回 数の節約を計ったコード（バンド行列法やスカイライン法）になっている． この計算結果は

ï

ï

ï

ï

ï

þ

ï

ï

ï

ï

ï

ý

ü

ï

ï

ï

ï

ï

î

ï

ï

ï

ï

ï

í

ì

EI

-=

ï

ï

ï

ï

ï

þ

ï

ï

ï

ï

ï

ý

ü

ï

ï

ï

ï

ï

î

ï

ï

ï

ï

ï

í

ì

b

b

b

b

11

65

.

32

75

.

10

75

.

21

5

.

9

5

.

11

25

.

6

416

.

3

8

p

w

w

w

w

5 5 4 4 3 3 2 2



である．ここで自由端の変位

w

₅は -32.65p/EI で, (2.17)において



=4 とおいて得 られる厳密解 -32p/EI に４要素でよく一致していることがわかる．

3.2 立体骨組解析

ここまでは簡単な直線梁について考えてきた．本節では任意の骨組構造の有限要素解 析を説明する．各梁要素が同一直線上になければ曲げ変形は伸び縮みや捩れになって伝わ っていく．立体梁要素の変位とそれに対応する力を図と表とに定義する．要素の座標系は x 軸に沿って梁を置くものとする．添え字の１は節点１の値であることを示す．

日本機械学会宇宙工学部門

u1 x方向の変位 X1 x方向の力 v1 y 方向の変位 Y1 y 方向の力 w1 z方向の変位 Z1 z方向の力 βx1 x 軸回りの捩れ角 Mx1 z 軸回りのモーメント βy1 y軸回りの回転 My1 y軸回りのモーメント βz1 z 軸回りの回転 Mz1 z 軸回りのモーメント

÷÷

ø

ö

çç

è

æ

¶

¶

=

b

¶

¶

=

b

¶

¶

=

b

x

w

,

x

w

,

y

w

z y x 図３．５ 立体梁要素

y

z

2

W

2

v

1

u

1

_u

2 ３次元梁要素を作るためには，x-z 平面における曲げ，x-y 平面における曲げ，x 軸方 向の伸び縮み，x 軸まわりの捩りを考える． まず，x-z 平面における曲げは式(3.5)より

ï

ï

þ

ï

ï

ý

ü

ï

ï

î

ï

ï

í

ì

b

w

b

w

ú

ú

ú

ú

û

ù

ê

ê

ê

ê

ë

é

-EI

=

ï

ï

þ

ï

ï

ý

ü

ï

ï

î

ï

ï

í

ì

2 y 2 1 y 1 2 2 2 2 3 2 y 2 1 y 1

L

4

L

6

L

2

L

6

L

6

12

L

6

12

L

2

L

6

L

4

L

6

L

6

12

L

6

12

L

M

Z

M

Z

(3.11) x-y 平面における曲げは

ï

ï

þ

ï

ï

ý

ü

ï

ï

î

ï

ï

í

ì

b

u

b

u

ú

ú

ú

ú

û

ù

ê

ê

ê

ê

ë

é

-EI

=

ï

ï

þ

ï

ï

ý

ü

ï

ï

î

ï

ï

í

ì

2 z 2 1 z 1 2 2 2 2 3 z 2 z 2 1 z 1

L

4

L

6

L

2

L

6

L

6

12

L

6

12

L

2

L

6

L

4

L

6

L

6

12

L

6

12

L

M

Y

M

Y

(3.12)

日本機械学会宇宙工学部門

下添字の y,z はそれぞれ y 軸，z 軸回りの値であることを示す．次に立体要素として１節 点６自由度を確保するために梁の伸縮要素（棒要素）と捩り要素をつくる．

3.2.1 棒要素と捩り要素

棒要素

図 3 ． 6  棒要素

EA

x

₁

x

₂

u

1

u

₂

L

(1.1)∼(1.3)より

( )

x

u

EA

E

A

A

F

¶

¶

=

e

=

s

=

0

=

¶

¶

x

u

の解を

u

=

c

₁

+

c

₂

x

として

u

(

0

)

=

u

₁、

u

(

L

)

=

u

₂の条件で 2 1

,

c

c

を求めると

x

L

)

u

u

(

u

u

=

₁

+

₂

-

₁ となる。よって

(

u

u

)

L

EA

F

=

₂

-

₁ で

x

₁

=

-

F

(

0

),

x

₂

=

F

(

L

)

とすると

þ

ý

ü

î

í

ì

ú

û

ù

ê

ë

é

-EA

=

þ

ý

ü

î

í

ì

2 1 2 1

1

1

1

1

u

u

L

X

X

(3.13) となる。

捩り要素

棒の場合と同様に(2.23)より

x

GJ

x

GJ

GJ

M

x T

_¶

b

¶

=

¶

q

¶

=

f

=

x

C

c

x

=

₁

+

₂

b

として

c

₁

,

c

₂を求める操作を行うと

GJ

図 3.7  棒の捩り要素

日本機械学会宇宙工学部門

þ

ý

ü

î

í

ì

b

b

ú

û

ù

ê

ë

é

-=

þ

ý

ü

î

í

ì

2 1 2 1

x

x

1

1

1

1

L

GJ

Mx

Mx

(3.14) を得る。

3.2.2 ３次元梁要素

式(3.11)(3.12)(3.13)(3.14)を総合して３次元梁要素の剛性方程式を組み上げると

[ ]

ï

ï

ï

ï

ï

ï

ï

ï

þ

ï

ï

ï

ï

ï

ï

ï

ï

ý

ü

ï

ï

ï

ï

ï

ï

ï

ï

î

ï

ï

ï

ï

ï

ï

ï

ï

í

ì

b

b

b

b

b

b

K

=

ï

ï

ï

ï

ï

ï

ï

ï

þ

ï

ï

ï

ï

ï

ï

ï

ï

ý

ü

ï

ï

ï

ï

ï

ï

ï

ï

î

ï

ï

ï

ï

ï

ï

ï

ï

í

ì

2 2 2 2 2 2 1 1 1 1 1 1 2 2 2 2 2 2 1 1 1 1 1 1

z

y

x

w

v

u

z

y

x

w

v

u

Mz

My

Mx

Z

Y

X

Mz

My

Mx

Z

Y

X

(3.15) ここに[K]は

ê

ê

ê

ê

ê

ê

ê

ê

ê

ê

ê

ê

ê

ê

ê

ê

ê

ë

é

-EI EI EI EI EI EI EI EI EA EI EI EI EI EI EI EI EI EA L 2 L 6 L 2 L 6 L GJ L 6 L 12 L 6 L 12 L L 4 L 6 L 4 L 6 L GJ L 6 L 12 L 6 L 12 L z 2 y 2 y 2 y 3 y 2 z 3 z z 2 y 2 2 y 3 y 2 z 3 z

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

日本機械学会宇宙工学部門

ú

ú

ú

ú

ú

ú

ú

ú

ú

ú

ú

ú

ú

ú

ú

ú

ú

û

ù

-EI EI EI EI EI EI EI EI EA EI EI EI EI EI EI EI EI EA L 4 L 6 L 4 L 6 L GJ L 6 L 12 L 6 L 12 L L 2 L 6 L 2 L 6 L GJ L 6 L L 6 L 12 L z 2 y 2 2 y 3 y 2 3 z z z y 2 y 2 y 3 y 2 z 3 z

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

12

0

0

0

0

0

0

0

0

0

0

0

(3.16) と与えられる．これが章末のプログラムの中のサブルーチン function[SK]=skmat(L,EA,EIZ,EIY,GJ)にプログラム化されている。 さて，多数の要素で構成される構造の解析においては，各要素が局所座標系で作られてい るのでその剛性行列を全体座標系に変換してから例３のような全体剛性行列を作る．変換 行列[T]により局所座標系と全体座標系（バー付きで表示）の関係は

[ ]

[ ]

[ ]

[ ]

ï

ï

ï

ï

ï

ï

ï

ï

þ

ï

ï

ï

ï

ï

ï

ï

ï

ý

ü

ï

ï

ï

ï

ï

ï

ï

ï

î

ï

ï

ï

ï

ï

ï

ï

ï

í

ì

b

b

b

b

b

b

ú

ú

ú

ú

û

ù

ê

ê

ê

ê

ë

é

T

T

T

T

=

ï

ï

ï

ï

ï

ï

ï

ï

þ

ï

ï

ï

ï

ï

ï

ï

ï

ý

ü

ï

ï

ï

ï

ï

ï

ï

ï

î

ï

ï

ï

ï

ï

ï

ï

ï

í

ì

b

b

b

b

b

b

z2 y2 x 2 z1 y x 1 1 e e e e z y x z y x

w

v

u

w

v

u

w

v

u

w

v

u

2 2 2 1 1 1 2 2 2 2 2 2 1 1 1 1 1 1 (3.17) ここに[

T

_e]は 3

´

3 掛けるの行列であり

日本機械学会宇宙工学部門

[ ]

ú

ú

ú

û

ù

ê

ê

ê

ë

é

-ú

ú

ú

û

ù

ê

ê

ê

ë

é

j

j

-j

j

=

T

q

mn/q

q

/

n

0

/q

m/q

n

m

cos

sin

0

sin

cos

0

0

0

1

e







(3.18) である．

Z

X

Y

x

z

y

図３．８ 局所座標と全体座標 ただし (

x

₁

,

y

₁

,

z

₁)を節点１の座標，(

x

₂

y

₂

z

₂)を節点２の座標として

(

x

₂

-

x

₁

)

/

L

m

=

(

y

₂

-

y

₁

)

/

L

=



(

)

2 1 2

z

/

L

q

m

z

n

=

-

=



+

である．部材の x 軸が全体座標系の z 軸と平行な場合には

[ ]

ú

ú

ú

û

ù

ê

ê

ê

ë

é

j

j

-j

j

=

T

0

cos

sin

0

sin

cos

1

0

0

e である．式(3.17)を

{ }

d

=

[ ]

T

{ }

d

(3.18) と書き表せば全体座標系での要素の剛性行列は式(3.26)より

[ ]

K

=

[ ] [ ][ ]

T

T

K

T

(3.20) となる．章末のプログラムではこの変換行列がサブルーチン function[h,fI]=henka 1 (xi,xj,yi,yj,zi,zj)で在る。式(3.20)の変換操作はメインプログラムの element

日本機械学会宇宙工学部門

assembly 部分で存される。 ここでの座標変換は実は情報不足で，実際にはＩ型梁，ｚ型梁のようにｙ方向曲げとｚ方 向曲げが異なるような場合の要素のｘ方向の回転を一意的に表示できない． ｘｙｚ系か らＸＹＺ系への変換には３軸の方向余弦，すなわち９個の情報が必要となるのにｉ，ｊ座 標の６個では不足である．これに対応するにはコードアングルの導入[3.4]や第３の点を指 定して向きを決める方式がある．

3.3 有限要素法解析の流れ

1. データ入力 (a) 節点数 (b) 要素数 (c) 節点の座標 (d) 要素を構成する節点の構成 (e) 各要素の EI,EA,GJ (f) 変位の境界条件 (g) 外力 2. 各要素での剛性行列の作成 3. 各要素の全体座標系への変換 4. 全体剛性行列に組み入れ 5. 変位境界条件の導入 6. 連立方程式を解く 7. 求められた変位から反力を求める 8. 求められた変位から要素の応力を求める

3.4 一般定式化

3.4.1 有限要素法における関係式

変位ー節点関数

{ } { }{ }

u

= N

d

(3.21) 梁の場合，式(3.3) 歪ー変位関係式

{ }

e

=

[ ]

B

{ }

d

(3.22)

日本機械学会宇宙工学部門

この関係式は変位節点関係式を微分することにより得られる．梁の場合 p.11 で説明した ように 2 2

dx

w

d

z

z

₌

-r

=

e

(

) (

) (

)

{

}

ï

ï

þ

ï

ï

ý

ü

ï

ï

î

ï

ï

í

ì

b

b

x

+

-x

-x

+

-x

+

-=

2 2 1 1 2

_w

w

3

1

L

6

3

3

2

L

6

3

L

2

z

応力ー歪関係式

{ }

σ =

[ ]

D

{ }

ε

(3.23) 梁要素の場合

s E

=

e

である．

3.4.2 仕事

歪エネルギー

{ } { }

e

s

=

{ }

s

[ ] [ ] [ ]

{ }

d

=

P

T T T

ò

ò

dV

B

D

B

dV

2

1

2

1

(3.24) ここに

[ ]

K

=

ò

[ ] [ ] [ ]

B

T

D

B

dv

(3.25) とすれば

{ }

d

[ ]

K

{ }

d

=

P

T

2

1

(3.26) 外部力による仕事

{ } { }

F

U

₂

=

d

T (3.27) 物体力による仕事

{ } { }

w

p

dV

{ } { } { }

N

p

dV

{ }

{ }

F

p

U

₁

=

ò

=

d

T

ò

T

=

d

T (3.28) よって物体力と等価な節点力

F

_pは

{ }

F

_p

=

-

ò

{ } { }

N

T

p

dV

(3.29)

日本機械学会宇宙工学部門

振動などによる慣性力も物体力で，時間に依存する力は

{ }

{ }

d

m

dt

u

d

dV

p

₂ 2

-=

(3.30) とおくことができる．したがって慣性力の場合

{ }

[ ]

{ }

2 2

dt

d

M

F

_p

=

d

(3.31) ここに

[ ]

M

=

ò

{ } { }

N

T

N

r

dV

(3.32) である．

3.4.3 エネルギー原理

物体力

p

と集中荷重

F

を受け，物体内部で変位

{ }

u

と節点変位

{ }

d

とを生じる構造物 の歪エネルギーを

P

とすると，荷重による外部仕事は歪エネルギーの変化に等しくなけれ ばならないので

( )

P

-

(

2

-

d

p

)

=

0

T

F

U

d

d

(3.33) すなわち与えられた変位に対して全ポテンシャルエネルギーが停留値をとる．歪エネルギ ー

P

は常に正の値をとるのでこの停留値は最小値である．これがポテンシャルエネルギー 最小の原理である． 系の平衡状態では変位に関してしてポテンシャルエネルギーは最小値をとるので，そのと きにはポテンシャルエネルギーを変位

d

Tで微分した値がゼロとなる．すなわち

[ ] { }

K

d

=

{ }

F

-

{ }

F

p (3.34) として剛性方程式が得られる．動的な運動を考えるときは

[ ]

K

{ }

d

( )

t

+

[ ]

M

{ }



d



( )

t

=

{

F

( )

t

}

(3.35) でこれは振動方程式となっており，先ほど定義した[M]は質量行列(mass matrix)と呼ば れる．

3.4.4 三角形平面要素

前節の有限要素法の定式化にしたがって演習問題として三角形平面要素を導いてみる． 1 2 3 y x

日本機械学会宇宙工学部門

変位−節点関係式

y

c

x

c

c

u

=

₁

+

₂

+

₃ (3.36a)

y

c

x

c

c

v

=

₄

+

₅

+

₆ (3.36b) とする．ここでは核節点で x 方向とｙ方向に２自由度、全３節点で計６個の自由度がある ので未定係数をｃ１∼ｃ６の６個としてある。これらは完全１次多項式となっている。点１ の(x,y)座標を(x1,y1)とし、そこの変位を u1,v1とすれば u₁ = c₁ + c₂x₁ +c₃y₁ u₂ = c₁ +c₂x₂ +c₃y₂ u₃ = c₁ + c₂x₃ +c₃y₃ なので

(

)

{

}

ï þ ï ý ü ï î ï í ì = 3 2 1 1 c c c y x y , x u であり、同様にv(x,y)も得られ逆行列を計算して

ú

û

ù

ê

ë

é

+

+

+

+

+

+

+

+

=

þ

ý

ü

î

í

ì

y

c

x

b

a

0

y

c

x

b

a

0

0

y

c

x

b

a

0

y

c

x

b

a

S

2

1

v

u

2 2 2 1 1 1 2 2 2 1 1 1

ï

ï

ï

ï

þ

ïï

ï

ï

ý

ü

ï

ï

ï

ï

î

ïï

ï

ï

í

ì

´

3 3 2 2 1 1

v

u

v

u

v

u

(3.37) を得る。ここに 2 3 3 2 1

x

y

x

y

a

=

-2 3 1 3 2 1

x

x

c

y

y

b

-=

-=

ï

þ

ï

ý

ü

ï

î

ï

í

ì

ú

ú

ú

û

ù

ê

ê

ê

ë

é

=

ï

þ

ï

ý

ü

ï

î

ï

í

ì

-3 2 1 1 3 3 2 2 1 1 3 2 1

1

1

1

u

u

u

y

x

y

x

y

x

c

c

c

日本機械学会宇宙工学部門

ú

ú

ú

û

ù

ê

ê

ê

ë

é

=

3 3 2 2 1 1

y

x

1

y

x

1

y

x

1

det

2

1

S

で S は三角形の面積となっている．

a

₁などの添え字はサイクリックに変わっていくもの とする． 歪−変位関係式は歪

{ }

e

が

{ }

ï

þ

ï

ý

ü

ï

î

ï

í

ì

+

=

ï

þ

ï

ý

ü

ï

î

ï

í

ì

g

e

e

=

e

¶ ¶ ¶ ¶ ¶ ¶ ¶ ¶ x u y u y v x u xy y x (3.38) として先ほどの変位関数を微分することにより

{ }

e

=

[ ]

B

{ }

d

ï

ï

ï

ï

þ

ïï

ï

ï

ý

ü

ï

ï

ï

ï

î

ïï

ï

ï

í

ì

ú

ú

ú

û

ù

ê

ê

ê

ë

é

=

ï

þ

ï

ý

ü

ï

î

ï

í

ì

g

e

e

3 3 2 1 1 1 3 3 2 2 1 1 3 2 1 3 2 1 xy y x

v

u

v

u

v

u

b

c

b

c

b

c

c

0

c

0

c

0

0

b

0

b

0

b

S

1

(3.39) が得られる． 薄い板が平面内で引張られるような平面応力状態では，応力−歪関係式は

{ }

s

=

[ ]

D

{ }

e

ï

þ

ï

ý

ü

ï

î

ï

í

ì

g

e

e

ú

ú

ú

û

ù

ê

ê

ê

ë

é

n

n

n

-E

=

ï

þ

ï

ý

ü

ï

î

ï

í

ì

t

s

s

n -xy y x xy y x 2 1 2

0

0

0

1

0

1

1

(3.40) である．よって剛性行列は

[ ] [ ] [ ][ ]

K

=

B

T

B

=

[ ] [ ][ ]

B

T

B

ò

D

dV

hS

D

(3.41) として積分なしで得られる．

日本機械学会宇宙工学部門

参考文献

有限要素法の具体的なプログラム作成について [3.1]三本木茂夫，吉村信敏：有限要素法による構造解析プログラム−考え方と解説−， コンピュータによる構造工学講座 I-1B, 培風館 1970. [3.2]小堀為雄，吉田博：有限要素法による構造解析プログラム，丸善,1980. [3.3]山田嘉昭監訳：CAL-90/SSTAN で解く有限要素構造解析（フロッピディスク付） 丸善， （12000 円） [3.4] 信原泰夫，桜井達美，吉村信敏：有限要素法のプログラムデザイン，コンピュータ による構造工学講座 Ⅱ-2B, 培風館 1970， pp.83-88.