73

有限要素法

•

概要

•

弱形式，要素分割，変位の有限要素基底による近似，ガラーキン 法の適用

•

仮想仕事の原理

•

弱形式と強形式

•

最小ポテンシャルエネルギーの原理

•

重み付き残差法と弱形式

•

ラグランジュ定数による変位規定境界条件の考慮

•

１次元問題の例

•

２次元問題の例

•

要素座標のパラメータと各種要素

•

ガウスの数値積分公式

•

応力の計算

有限要素法

•

概要

•

弱形式，要素分割，変位の有限要素基底による近似，ガラーキン 法の適用

•

仮想仕事の原理

•

弱形式と強形式

•

重み付き残差法と弱形式

•

１次元問題の例

•

２次元問題の例

•

要素座標のパラメータと各種要素

•

ガウスの数値積分公式

•

応力の計算

有限要素法の流れ

弱形式 仮想仕事の原理，最小ポテンシャ ルエネルギーの原理，重み付き残 差法，ガラーキン法 要素分割 （要素形状） 形状関数（有限要素基底） 連立方程式 解 ! ! (∇u · ∇w + bw) d! = ! "_t ¯tw d" u = N1u1 _{+ N}2u2 _{+ N}3u3 _{+ N}4u4 _{+ · · ·} [K ] {u} = {F}

仮想仕事の原理

釣合状態（平衡状態）において，幾何学的境界条件（変位規程境界条件を満足 する変位 u （可容変位）の変分δu（仮想変位）とすると，釣合状態にあり，力 学的境界条件を満足している応力（可容応力）が，仮想変位に対してなす仕事 （内部仕事）は，外力と物体力が仮想変位に対してなす仕事（外部仕事）に等 しい．このことを仮想仕事の原理という． x L dx σ _{σ + dσ} δu bδu + d(δu) 左端が固定された断面積A，長さLの棒の外力が仮想変位に対してなす仕事はW は次のようになる． 左端では，仮想変位δu = 0であるから，左端の外力の仕事 (= 0)を追加して変形 すると次のようになる． W = ! _L 0 b δu dV + Aσ(L)δu(L) W = A ! _L

0 b δu dx + Aσ(L)δu(L) − Aσ(0)δu(0)

σ (L)

（式変形の続き）          釣合方程式より= 0 積分の形に戻す W = A ! _L

0 b δu dx + Aσ (L)δu(L) − Aσ(0)δu(0) = A

! _L 0 b δu dx + A ! _L 0 d(σ δu) dx dx = A ! _L 0 b δu dx + A ! _L 0 dσ dx δu dx + A ! _L 0 σ d(δu) dx dx = A ! _L 0 "_dσ dx + b # δu dx + A ! _L 0 σ d(δu) dx dx = A ! _L 0 σ d(δu) dx dx = A ! _L 0 σ δε dx δuに対応するひずみ 外力が仮想変位に対してなす仕事 は内部応力がなす仕事に等しい A ! _L σ δε dx = A ! _L b δu dx + Aσ(L)δu(L) 仮想仕事の原理

仮想仕事の原理（一般の場合）

W = ! ! (b_xδu + b_yδv + b_zδw) d! + ! #_t (t_xδu + t_yδv + t_zδw)d# = ! ! (b_xδu + b_yδv + b_zδw) d! + ! #_t+#_u (t_xδu + t_yδv + t_zδw) d# = ! ! (b_xδu + b_yδv + b_zδw) d! + ! #{(σ xxnx + σyxny + σzxnz)δu + (σxynx + σyyny + σzynz)δv + (σxznx + σyzny + σzznz)δw} d# = ! ! (b_xδu + b_yδv + b_zδw) d! + ! #{(σ xxδu + σxyδv + σxzδw)nx + (σyxδu + σyyδv + σyzδw)ny + (σzxδu + σzyδv + σzzδw)nz} d# （次ページに続く） 物体力の仮想仕事 表面力の仮想仕事 変位規定の境界も追加 コーシーの公式により表面力を応力で表現 外向き単位法線ベクトルの成分でまとめる

（式変形の続き） W = · · · = ! ! (b_xδu + b_yδv + b_zδw)d! + ! ! " ∂(σ_xxδu + σ_xyδv + σ_xzδw) ∂x + ∂(σ_yxδu + σ_yyδv + σ_yzδw) ∂y + ∂(σzxδu + σzyδv + σzzδw) ∂z # d! = ! ! (b_xδu + b_yδv + b_zδw)d! + ! ! "$ ∂σ_xx ∂x + ∂σ_yx ∂y + ∂σ_zx ∂z % δu + $ ∂σ_xy ∂x + ∂σ_yy ∂y + ∂σ_zy ∂z % δv + $ ∂σ_xz ∂x + ∂σ_yz ∂y + ∂σ_zz ∂z % δw # d! + ! ! " σ_xx ∂(δu) ∂x + σyy ∂(δv) ∂y + σzz ∂(δw) ∂z + σxy $ ∂(δu) ∂y + ∂(δv) ∂x % + σyz $ ∂(δw) ∂y + ∂(δv) ∂z % + σzx $ ∂(δw) ∂x + ∂(δu) ∂z %# d! 発散定理により領域積分に変換 法線ベクトルはその方向の偏微分にし， 境界積分を領域積分に変える． 仮想変位の成分でまとめる 残りの仮想変位の微分の項をまとめる まとめる際に応力テンソルの対称性を利用 垂直ひずみ

（式変形の続き） 剪断応力，工学剪断ひずみの表現を利用 釣合い方程式よりこれらは０ したがって，一般の場合の仮想仕事の原理は次のようになる． これを簡潔に書くと次のようになる． W = · · · = ! ! "# ∂σ_xx ∂x + ∂σ_yx ∂y + ∂σ_zx ∂z + bx $ δu + # ∂σ_xy ∂x + ∂σ_yy ∂y + ∂σ_zy ∂z + by $ δv + # ∂σ_xz ∂x + ∂σ_yz ∂y + ∂σ_zz ∂z + bz $ δw % d! + ! ! (σ_xδε_x + σ_yδε_y + σ_zδε_z + τ_xyδγ_xy + τ_xyδγ_xy + τ_xyδγ_xy) d! = ! ! (σ_xδε_x + σ_yδε_y + σ_zδε_z + τ_xyδγ_xy + τ_xyδγ_xy + τ_xyδγ_xy) d! ! ! (σ_xδε_x + σ_yδε_y + σ_zδε_z + τ_xyδγ_xy + τ_xyδγ_xy + τ_xyδγ_xy) d! = ! ! (b_xδu + b_yδv + b_zδw) d! + ! '_t (t_xδu + t_yδv + t_zδw) d' ! !!δε"{σ } d! = ! !!δu"{b} d! + ! %_t!δu"{t} d%

仮想仕事の原理（弱形式） ! ! (σ_xδε_x + σ_yδε_y + σ_zδε_z + τ_xyδγ_xy + τ_xyδγ_xy + τ_xyδγ_xy) d! = ! ! (b_xδu + b_yδv + b_zδw) d! + ! '_t (t_xδu + t_yδv + t_zδw)d'

弱形式と強形式

仮想仕事の原理の導出過程で，釣合い方程式および力学的境界条件が用いら れている．また，仮想仕事の原理から部分積分の公式を逆向きに適用して変 形していくと，釣り合い方程式と力学的境界条件が得られる．したがって， 仮想仕事の原理は，釣り合い方程式と力学的境界条件に等価である． 仮想仕事の原理の式を一般に弱形式という． または ! !!δε"{σ } d! = ! !!δu"{b} d! + ! %_t!δu"{t} d%

最小ポテンシャルエネルギーの原理

応力成分がひずみエネルギー密度関数の導関数として表されるものとする. 等方線形弾性体の場合のひずみエネルギー密度関数は，次のようになる． たとえばU0のεxに対する導関数がσxとなっていることが次のように確かめ られる． また，表面力や物体力が保存力であり，それらのポテンシャルをVT，VBと すると， {σ } = ! ∂U₀ ∂ε " U_{0 =} 1 2(σxεx + σyεy + σzεz + τxyγxy + τyzγyz + τzxγzx) = _{2(1−2ν)(1+ν)}νE (ε_x + ε_y + ε_z)2 + _2(1+ν)E ! ε_x2 + ε2_y + ε2_z + 1₂γ_xy2 + 1₂γ_yz2 + 1₂γ_zx2 " ∂U₀ ∂ε_x = νE (1 − 2ν)(1 + ν)(εx + εy + εz) + E 1 + ν εx = σx {t} = ! −∂VT ∂t " , {b} = ! −∂VB ∂b "

したがって，仮想仕事の原理は次のようになる． ここで， ! !!δε" " ∂U₀ ∂ε # d! = ! !!δu"{− ∂V_B ∂u } d! + ! %_t!δu"{− ∂V_T ∂u } d% ! ! δU₀ d! = − ! ! δV_B d! − ! #_t δV_T d# ! = ! " U_{0 d" +} ! " V_{B d" +} ! #_t V_T d# δ !" " U_{0 d" +} " " VB d" + " #_t VT d# # = 0 を全ポテンシャルエネルギーという．物体力と表面力が死荷重（仮想変位 を与えても値が変化しない荷重）のときは，VBとVTはそれぞれ となり，Πは次のようになる． ! = ! " U_{0 d" −} ! ""u#{b} d" − ! #_t"u#{t} d# = ! U _{d" −} ! b u _{d" −} ! t u d# VB = −"u#{b} = −biui, VT = −"u#{t} = −tiui U0の変分

最小ポテンシャルエネルギーの原理

Ω内で連続で，Γ_u上で変位の境界条件を満足するすべての許容される変位

成分の中で，_{正解となるものは全ポテンシャルエネルギーΠを最小にす}

る_{．詳細は省略するが，δΠ = 0 とすると，}

!" = "(u + δu) − "(u) = δ" + δ2" = δ2" > 0

が示せて， !(u + δu) > !(u) となることが分かる． 特に線形弾性体の場合， U_{0 =} 1 2 !!"{σ } = 1 2 !ε" [D] {ε} であるから， ! = 1 2 ! "!ε" [D] {ε} d" − ! "!u"{b} d" − ! $_t!u"{t} d$ よって，δΠ = 0 は次のように仮想仕事の原理と同じ弱形式が得られる． δ" = ! #!δε" [D] {ε} d# − ! #!δu"{b} d# − ! %_t!δu"{t} d% = ! #!δε"{σ } d# − ! #!δu"{b} d# − ! %_t !δu"{t} d% = 0

重み付き残差式と弱形式

微分方程式を直接解く代わりに，近似解を微分方程式に代入したときの誤 差に重みをかけ，領域全体で積分したものを０となるように近似解を構成 する方法を，重み付き残差法という．重み付き残差法の式を出発点とし て，部分積分を行うと有限要素法で用いる弱形式が得られる． 以下，簡単のため添字記号を用いて説明する． Γuで0となる を重み関数とするとし，次の積分を考える．wj 左辺第１項を部分積分すると次にようになる． ! ! t_ju∗_{j d! −} ! " σ_{i j}u∗_j,i d" = ! ! t_ju∗_{j d! −} 1 2 ! " σ_{i j}(u∗_j,i + u∗_j,i) d" = ! ! t_ju∗_{j d! −} 1 2 ! " (σ_{i j}u∗_j,i + σ_jiu∗_j,i) d" = ! ! t_ju∗_{j d! −} 1 2 ! " (σ_{i j}u∗_j,i + σ_{i j}u∗_{i, j}) d" = ! !_u+!_t t_ju∗_{j d! −} ! " σ_{i j}ε_{i j}∗ d" = ! !_t t_iu∗_{i d! −} ! " σ_{i j}ε_{i j}∗ d" ! ! (σ_{i j,i} + b_j)u∗_j d! = 0

したがって，下の式は次のようになる． ! !_t t_iu_i∗ _{d! −} ! " σ_{i j}ε_{i j}∗ d" + ! " b_ju∗_{j d" = 0} この式の左辺第１項のt iは，近似解に関係づけられる表面力であるので， この部分を力学的境界条件に一致するようすするために，重み付き残差式 として，次の形を出発点とする． ! ! (σ_{i j,i} + b_j)u∗_j d! − ! #_t (t_j − ¯t_j) d# = 0 この式は，領域内部で近似解による微分方程式の残差と境界上の残差を０ とする形になっている．この式の左辺を部分積分すると，次の弱形式が得 られる． 重み関数として，仮想変位を考えれば，仮想仕事の原理から導出される弱 形式と一致していることが分かる． ! ! σ_{i j}ε_{i j}∗ d! = ! ! b_iu_i∗ _{d! +} ! $_t ¯tiu_i∗ d$

ラグランジュ乗数を用いた変位規定境界条件の組み込み

最小ポテンシャルエネルギーの原理において，ラグランジュ乗数を用いて 変位規定境界条件を組み込む．すなわち， ! = 1 2 ! " σ_{i j}ε_{i j} d" − ! " b_iu_{i d" −} ! %_t ¯tiui d% − ! %_u λ_i(u_i − ¯u_i) d% 第１変分をとると， δ" = 1 2 ! # σ_{i j}δε_{i j} d# − ! # b_iδu_i d# − ! &_t ¯tiδui d& − ! &_u δλ_i(u_i − ¯u_i) d& − ! &_u λ_iδu_i d& ここで，部分積分により， ! ! σ_{i j}δε_{i j} d! ! %_u+%_t t_iδu_i d% − ! ! σ_{i j,i}δu _j d! よって， δ" = − ! # (σ_{i j,i} + b_j)δu _j d# − ! # b_iδu_i d# + ! %_t (t_i − ¯t_i)δu_i d% − ! %_u (λ_i − t_i)δu_i d% − ! %_u δλ_i(u_i − ¯u_i) d% Γ_u で δu_i = 0 でなくてよい．

（続き） δΠ=0の条件より， σ_{i j,i} + b_j = 0 t_{i = ¯ti} λ_i = t_i, u_i = ¯u_i ! （ 内で） !_t （  上で） !_u （  上で） すなわち，ラグランジュ定数はΓu上で表面力の意味を持ち，変位規定境界 条件を満足するようにできることが分かる． δΠの第１項を部分積分して得られる弱形式は，次のようになる． δ" = − ! # (σ_{i j,i} + b_j)δu _j d# + ! %_t (t_i − ¯t_i)δu_i d% − ! %_u (λ_i − t_i)δu_i d% − ! %_u δλ_i(u_i − ¯u_i) d% ! ! σ_{i j}δε_{i j d!} = ! ! b_jδu _j d! + ! %_t ¯tiδui d% + ! %_u (λ_i − t_i)δu_i d% + ! %_u δλ_i(u_i − ¯u_i) d%

弱形式に基づく有限要素法

（１次元問題の例）

強形式 x L σ (L) 断面積A 重み付き残差式 dσ dx + b = 0, u(0) = 0, σ(L) = F/A ! _L 0 w " dσ dx + b # Adx = 0 １回部分積分することにより弱形式が得られる． Hookeの法則，変位・ひずみ関係式を用いると次のようになる． ! _{L dw du " #L} ! _L ! wσ A "L 0 − # _L 0 dw dx σ Adx + # _L 0 b Adx = 0

有限要素法

方針 領域を要素に分割し，積分を要素ごとに評価する （この例では区間に分割） 解を形状関数を用いて節点値の線形結合として近似 重み関数（仮想変位）に同じ基底関数を用いる 連立方程式 弱形式を求める ! _L 0 dw dx E du dx Adx = " wσ A#L 0 + ! _L 0 b Adx = 0 有限要素基底の導入 ガラーキン法

要素分割と変位の近似

次のように，区間をいくつかの小区間（_{ここでは幅Δで4等分割}）に分割し て，積分を区間ごとに評価する． = L u 0 Δ 2Δ 3Δ 4Δ x u1 u2 u3 u4 u5 uを下の図のように区間ごとにuの近似値（節点値）を直線で結んだものと して近似する． ! e " !_e dw dx E du

dx Adx = w(L)σ(L)A − w(0)σ(0)A +

!

e

"

このとき，uは次のように書くことができる． u(x) ≈ 5 ! a=1 Na(x)ua Na(x) は形状関数と呼ばれ，この近似の場合次のような関数である． 0 x ua 1 xa xa+1 xa–1 Na(x) 形状関数は，考えている節点で1，考えている節点から隣接する節点まで連 続的に，0に変化し，考えている節点を含まない要素ではすべての点で0と なる関数である．

要素 e に着目すると，u(x) は要素 e 内のでは次のように書ける．

u = Ne(x)ue + Ne+1(x)ue+1

要素 e 内では，Ne_(x)とNe+1_{(x) は次の図のようになっている．} 0 x 1 xe+1 xe Δ Ne+1(x) Ne(x) このとき要素 e 内でのu(x)は下の図のように変化する． Ne+1(x) = x − x e xe+1 _{− x}e = x − xe ! Ne(x) = 1 − x − x e xe+1 _{− x}e = xe+1 _{− x} ! 0 1 xe+1 xe Δ Ne+1(x) Ne(x) u ue(x) ue+1(x) x

弱形式の重み関数を考える． 重み関数には，ある特定の節点の仮想的な変位（仮想変位）を用いる．こ の節点を節点aとすると， このとき，節点aを含まない要素の積分はすべて0となり，節点aを含む要 素のみ積分を計算すればよいことが分かる．  逆に，ある要素 _e の積分を行う場合，要素内の節点に仮想的変位を与え る場合のみが積分の計算対象となることがわかる．  したがって，この例の場合は要素eの積分に対しては仮想変位を与えるの は節点eの場合と節点e+1の場合ということになる． 要素Δeの積分は次のようになる． w = w(x) = Na(x)wa すなわち， ! e " !_e dw dx E du

dx Adx = w(L)σ(L)A − w(0)σ(0)A +

! e " !_e b Adx = 0 wa ! e " _xe_+! xe d Na dx E du

dx Adx = waNa(L)σ (L)A − waNa(0)σ (0)A + wa

! e " _xe_+! xe N a_{b Adx = 0} ! e " _xe_+! xe d Na dx E du dx dx = Na(L)σ (L) − Na(0)σ (0) + ! e " _xe_+! xe N a_{b dx = 0}

要素e内の変位の導関数は次のようになる． du dx = d Ne dx u e + d N e+1 dx ue+1 = ! d Ne dx d Ne+1 dx " ! ue ue+1 " ここで d Ne dx = −1 ! , d Ne+1 dx = 1 ! であるから， まず      のときを考えるとw = Newe

ke,e ₌ E , ke,e = −E

! _xe_+! xe d Ne dx E du dx dx = ! _xe_+! xe d Ne dx E " d Ne dx d Ne+1 dx # $ ue ue+1 % dx = "! _xe_+! xe d Ne dx E d Ne dx dx ! _xe_+! xe d Ne dx E d Ne+1 dx dx # $ ue ue+1 %

≡ &k₁e,e k₂e,e' $

ue ue+1

ただし， 物体力の積分についても同様に考えると，     のとき 同様に      _{w = N}e+1_we+1  _{のときは次のようになる．} ! _xe_+! xe dw dx E du dx dx = " k₁e+1,e k₂e+1,e# $ ue ue+1 % k₁e+1,e ₌ −E ! , k e+1,e 2 = E ! ! _xe_+! xe N e_{b dx =} "! 3 ! 6 # $ be be+1 % = Se,e w = Newe        のとき_{w = N}e+1_we+1 ! _xe_+! xe N e_{b dx =} "! 6 ! 3 # $ be be+1 % = Se+1,e

以上をまとめると，仮想変位 we _{(e=2, 3, 4) に対しては，} すなわち， 左端の要素を積分するときには，積分結果が 0 にならないのは w = N1_w1 _ま たは w = N2_w2 _{の場合のみである．w = N}1_w1 _{のときは x = 0 で N}1_{(0) = 1 とな} るから，右辺の (0)の項は考える必要があるが，N5_{(L)=0であるから (L)} の項は0となる． 右端 x = L を含む要素（この場合はΔ4）では，w = N4_w4 または w = N5_w5 の 場合に積分結果が 0 とならない．特に w = N5_w5 _{のときは， x = L で N}5 _{= 1} となるから，右辺の (L) の項を考える必要がある． w = N2w2, N3w3, N4w4 のときは，右辺の応力に関する項は 0 となる． σ σ σ ! k₁e,e−1 k₂e,e−1" # ue−1 ue $

+ %k₁e,e k₂e,e& # ue ue+1 $ − (Se,e−1 + Se,e) = 0 !

k₁e,e−1 k₂e,e−1 _{+ k1}e,e k₂e,e"

     ue−1 ue ue+1      = S e,e−1 _{+ S}e,e

左端 x = 0に仮想変位を与えるとき，すなわちw = N1_w1 _{に対しては，次のよ} うになる． 以上をまとめると， ! k₁5,4 k₂5,4" # u4 u5 $ = S5,4 + σ (L) 左端 x = 0に仮想変位を与えるとき，すなわちw = N5_w5 _{に対しては，次のよ} うになる．          E ! −E! 0 0 0 −E ! 2E ! −E! 0 0 0 −E_! 2E_! −E_! 0 0 0 −E_! 2E_! −E_! 0 0 0 −E_{! !}E                         ¯ u₁ u₂ u₃ u₄ u₅                =                S1,1 _{− σ (0)} S2,1 _{+ S}2,2 S3,2 _{+ S}3,3 S4,3 _{+ S}4,4 S5,4 _{+ ¯σ (L)}                ! k₁1,1 k₂1,1" # u1 u2 $ = S1,1 − σ (0)

   は左端の変位であり，変位規定境界条件より既知量である．したがっ て，(対応する列) ×  を右辺に移項するとともに，1行目を別に分けて書 くと次のようになる． ¯u1 ¯u1       2E ! −E! 0 0 −E ! 2E ! −E! 0 0 −E_! 2E_! −E_! 0 0 −E_{! !}E                u₂ u₃ u₄ u₅          =            S2,1 _{+ S}2,2 _{+ !}E _¯u1 S3,2 _{+ S}3,3 S4,3 _{+ S}4,4 S5,4 _{+ ¯σ (L)}            すなわち [K ] {u} = {F} および この式は，すべての節点の変位が得られたあとで，変位既知の左端の応力 ! _E ! −E! 0 0 0 "              ¯u1 u₂ u₃ u₄ u₅              = S1,1 − σ (0)

100 これを解くと，正解が次のように得られる． 特に，_{¯u1 = 0}，物体力が存在しないときは，次のようになる． ! ! wq_i,i d! _{= 0} ! " wq_in_i d" − ! ! w_iq_i d! = 0 ! ! w_iq_i d! = ! " wq_in_i d" q_{i = −ki} ∂T ∂x_i ! ! w_,i " −ki ∂T ∂x_i # d! ₌ ! " wq_n d"

For two-dimensional case, ! ! $ ∂w ∂x " −kx ∂T ∂x # + ∂w ∂y " −ky ∂T ∂y #% d! ₌ ! " wq_n d" ! ! " ∂w ∂x ∂w ∂y # & −kx 0 0 _−ky '          ∂T ∂x ∂T ∂y          d! ₌ ! " wq_n d" T ₌ 3 / b₌₁ NbT b ∂T ∂x = 3 / b₌₁ ∂Nb ∂x T b ! $ w_i " −ki ∂T ∂x_i # d$ ₌ ! " wq_in_i d" 7 E !      2 ₋₁ 0 0 −1 2 −1 0 0 ₋₁ 2 ₋₁ 0 0 ₋₁ 1               u₂ u₃ u₄ u₅          =          0 0 0 ¯σ (L)          u_{2 =} ¯σ (L)" E , u3 = 2 ¯σ(L)" E , u4 = 3 ¯σ(L)" E , u5 = 4 ¯σ(L)" E また，x=0の応力（反力）は次のように得られる． σ (0) = − ! _"E −E_" 0 0 0 "              ¯u1 u₂ u₃ u₄ u₅              = ¯σ (L)

変位規定境界条件を弱形式に組み込む場合

! _L 0 σ δε Adx = ! _L 0 bδu Adx + " Aσ (x)δu(x)#L 0 ラグランジュ乗数を使って，変位規定境界条件を制約条件として組み込む 方法は，結局次の弱形式で変位規定境界条件が科されている境界において も0でない仮想変位を考えるかわりに，その境界の節点値に境界条件値を 規定することに帰着する． この式から次の連立方程式が得られる． 右辺に移項 x = 0 の節点反力を求めるために使用 AE !          1 ₋₁ 0 0 0 −1 2 −1 0 0 0 ₋₁ 2 ₋₁ 0 0 0 ₋₁ 2 ₋₁ 0 0 0 ₋₁ 1                           ¯u1 u₂ u₃ u₄ u₅                  =                  −Aσ (0) 0 0 0 Aσ (L)                 

結局，次の2組の式に帰着する． AE !       2 ₋₁ 0 0 −1 2 −1 0 0 ₋₁ 2 ₋₁ 0 0 ₋₁ 1                  u₂ u₃ u₄ u₅            =            AE ¯u1 ! 0 0 Aσ (L)            −Aσ (0) = AE " "1 − 1 0 0 0#            u₁ u₂ u₃ u₄ u₅            (i) 未知変位を求める式 (ii) 変位規定境界における節点反力を求める式 (a) (b) 未知変位を式(a)より求めた後，それらを式(b)に代入すると，x=0の応力が 次のように得られる． σ (0) = σ(L)

まとめ

•

要素ごとに積分

•

要素ごとの積分結果が0とならない場合は，仮想変位

（重み関数）がその要素の節点変位の時のみ．積分の

結果は，全体のマトリックスの，仮想変位の節点番号

に対応する行，積分を計算する要素の節点番号に対応

する列に足し込む．

•

物体力などから計算される等価節点力は，右辺のベク

トルの仮想変位に対応する行に足し込む．

•

表面力が規定されている境界の節点を含む要素を積分

するときは，仮想変位が境界節点に対応するときは境

界積分を計算し，得られた等価節点力を対応する行に

足し込む．

２次元弾性問題における有限要素法

弱形式 ただし， [D] = E (1 + ν)(1 − 2ν)   1 − ν ν 0 ν 1 − ν 0 0 0 1−2ν₂      ε_x ε_y γ_xy    =      ∂ ∂x 0 0 _∂∂_y ∂ ∂y ∂∂x      -u v . ! ! (σ_xδε_x + σ_yδε_y + τ_xyδγ_xy) d! = ! ! (b_xδu + b_yδv) d! + ! '_t (t_xδu + t_yδv) d' σ_xδε_x + σ_yδε_y + τ_xyδγ_xy = !δε_x δε_y δγ_xy" [D]    ε_x ε_y γ_xy   

要素分割（三角形）

弱形式 形状関数を用いた変位の近似 ! u v " = N1 ! u1 v1 " + N2 ! u2 v2 " + N3 ! u3 v3 " + · · · Na a 1 ! e " !e!δε x δεy δγxy" [D]    ε_x ε_y γ_xy    d% = ! e " !e!δu δv" ) b_x b_y * d% + ! k " &_tk!δu δv" ) t_x t_y * d&

それぞれの三角形に対しては，形状関数はそれぞれ次のような分布そして いる． 1 2 3 1 1 3 2 1 1 2 3 1

三角形１次要素の形状関数の具体形

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

要素ごとの積分の計算

形状関数の導関数 ∂ N1 ∂x = B₁ 2", ∂ N1 ∂y = C₁ 2" ∂ N2 ∂x = B₂ 2", ∂ N2 ∂y = C₂ 2" ∂ N3 ∂x = B₃ 2", ∂ N3 ∂y = C₃ 2"    ε_x ε_y γ_xy    =        ∂N1 ∂x 0 0 ∂ N_∂y1 ∂N1 ∂y ∂N 1 ∂x        -u1 v1 . +        ∂ N2 ∂x 0 0 ∂ N_∂y2 ∂ N2 ∂y ∂ N 2 ∂x        -u2 v2 . +        ∂ N3 ∂x 0 0 ∂ N_∂y3 ∂ N3 ∂y ∂ N 3 ∂x        -u3 v3 . = 3 / j=1        ∂ N j ∂x 0 0 ∂ N_∂yj ∂ N j ∂y ∂ N j ∂x        -u j v j . = 3 / j=1 0 B j1 -u j v j .























一定

要素ごとの積分の計算（続き）

! !e!δε x δεy δγxy" [D]    ε_x ε_y γ_xy    d% = !δu i _δvi_" (3 j=1 )! !e * Bi+*D+*B j+ d% , -u j v j . = !δui δvi" 3 (_* ki j+ -u j v j . !δεx δεy δγxy" =! δui δ vi"!Bi"T    ε_x ε_y γ_xy    = 3 ' j=1 ( B j) * u j v j +

+

三角形線形要素 のときは一定値 ひずみ，仮想ひずみを形状関数で表すと，結局次のようになる． ! e " !e!δεx δε_y δγ_xy" [D]    ε_x ε_y γ_xy    d% = ! e " !e!δu δv" ) b_x b_y * d% + ! k " &_tk!δu δv" ) t_x t_y * d&

領域積分が０でない値となるとき

e k i j e k i j e k i j

左図の赤丸ように，仮想変位が積分を評

価する要素外の節点の変位のときは，要

素 eの領域積分の値は0となる．

（仮想変

位の形状関数の値が要素 e内では常に0と

なるから．）

e k i j

仮想変位を与える節点（上図の赤丸）が積分を評価す

る要素内の節点のいずれかの節点の変位のとき

要素の局所節点番号と大域的節点番号

要素の局所節点番号には，対応する大域的節点番号がある．

要素ごとの積分の計算により，要素ごとの局所節点数の積分

結果が得られる．これをマトリックスの大域的節点番号に対

1 2 3 e k i j e

要素 e の局所節点番号

要素 e の大域的節点番号

表面力の境界積分の評価

仮想変位を与える節点が境界上の節点と一致すると きは，その節点を含む境界の積分も計算しなければ ならない． 線分境界上の点の座標は，パラメータを用いて次のように書ける． 図の点 i で仮想変位を考え，i j 間の線分境界上の境界 積分を評価する．この辺をΔStとすると x = xi + 1 + ξ₂ (xj − xi), −1 ≤ ξ ≤ 1 表面力も同様に直線上に変化するものとすると， t = ti + 1 + ξ₂ (tj − ti), −1 ≤ ξ ≤ 1 また，積分をパラメータξで評価するものとすると， d! = !"_dx dξ #2 + " dy_dξ #2 dξ = L₂ dξ i j k L ! !"_t (δu t_x + δv t_y) d" = ! !"_t!N i_δui _Ni_δvi_" " tx ty # d"

Ni(x, y) は境界上では節点 i で1，節点 j で 0 となり，節点 i と j を直線で結 んだ形となるので，パラメータξを用いると Ni(x, y) = 1 − ξ 2 となる．以上より， ! !"_t (δu t_x + δv t_y) d" = !δui δvi" ! ₁ −1      & 1−ξ 2 '2 t_xi ₊ 1−ξ₄ 2t_xj & 1−ξ 2 '2 t_yi ₊ 1−ξ₄ 2t_yj      L 2 dξ = !δui δvi"    L 3 txi + L6 txj L 3 tyi + L6 tyj    ! " !"_tk (δu t_x + δv t_y) d" = !δui δvi"      ' ( _L 3 txi + L6 txj ) ' ( _L 3 tyi + L6 tyj )      = !δu i _δvi_" -F_xi F_yi .

全体の連立方程式の組み立て

離散化された弱形式 ! e " !e!δε x δεy δγxy" [D]    ε_x ε_y γ_xy    d% = ! " !&_t !δu δv" ) tx ty * d& ! " !"_tk (δu t_x + δv t_y) d" = !δui δvi" # F_xi F_yi $ !δu δv"   n # j=1 $ ki j% & u j v j ' − ( F_xi F_yi ) = 0   ! e " !e!δε x δεy δγxy" [D]    ε_x ε_y γ_xy    d% = !δu i _δvi_" !n j=1 ) ki j* + u j v j , K {u} = {F} 積分の評価 代数方程式

形状関数のパラメータ表示

x y x1 x2 x3 x4 0 −1 1 1 −1 x1(−1, −1) x2(1, −1) x3(1, 1) x4(−1, 1) ξ η ξ u2 u3 u4 u1 0 −1 1 1 −1 η

∂Nb ∂ξ = ∂ Nb ∂x ∂x ∂ξ + ∂ Nb ∂y ∂y ∂ξ ∂Nb ∂η = ∂ Nb ∂x ∂x ∂η + ∂ Nb ∂y ∂y ∂η      ∂ Nb ∂ξ ∂ Nb ∂η      = J      ∂ Nb ∂x ∂ Nb ∂y      ２次元の要素の形状関数の実座標による導関数とパラメータによる導関数の 関係を求める．      ∂ Nb ∂x ∂ Nb ∂y      = J −1      ∂ Nb ∂ξ ∂ Nb ∂η      ! ka,b_j " ₌ # !"b ! Ba"T!D"!B j" d" d! = ! ! ! ! ∂x ∂ξ × ∂x ∂η ! ! ! ! dξdη 積分をパラメータで実行する ための面積要素







∂Nb ∂ξ ∂Nb ∂η







=





∂x ∂ξ ∂y ∂ξ ∂x ∂η ∂y ∂η











∂Nb ∂x ∂Nb ∂y







3次元のソリッド要素

ξ η ζ 1 2 3 4 5 6 7 8 N1 ₌ 1 8(1 − ξ)(1 − η)(1 − ζ), N2 = 1 8(1 + ξ)(1 − η)(1 − ζ) N3 ₌ 1 8(1 + ξ)(1 + η)(1 − ζ), N4 = 1 8(1 − ξ)(1 + η)(1 − ζ) N1 ₌ 1 8(1 − ξ)(1 − η)(1 + ζ), N2 = 1 8(1 + ξ)(1 − η)(1 + ζ) N3 ₌ 1(1 + ξ)(1 + η)(1 + ζ), N4 = 1(1 − ξ)(1 + η)(1 + ζ) ８節点六面体要素の形状関数

             ∂Nb ∂x ∂Nb ∂y ∂Nb ∂z              = J−1              ∂Nb ∂ξ ∂Nb ∂η ∂Nb ∂ζ              d! = ! ∂x ∂ξ × ∂x ∂ξ " · ∂x ∂η dξdηdζ = det J dξdηdζ 積分をパラメータで実行するための体積要素 形状関数の導関数              ∂ Nb ∂ξ ∂ Nb ∂η ∂ Nb ∂ζ              =        ∂x ∂ξ ∂y ∂ξ ∂ξ∂z ∂x ∂η ∂y ∂η ∂η∂z ∂x ∂ζ ∂y ∂ζ ∂ζ∂z                     ∂ Nb ∂x ∂ Nb ∂y ∂ Nb ∂z              = J              ∂ Nb ∂x ∂ Nb ∂y ∂ Nb ∂z             

数値積分（ガウスの数値積分公式）

! ₁ −1 f (ξ) dξ = n " i=1 f (ξ_i)w_i ξ_i, w_i (i = 1, 2, · · · , n) はそれぞれ積分点の座標と重み． n 0 1 2.000 000 000 000 000 0.577 350 269 189 626 2 1.000 000 000 000 000 0.774 596 669 241 483 0.000 000 000 000 000 3 0.555 555 555 555 556 0.888 888 888 888 889 0.861 136 311 594 053 0.339 981 043 584 856 4 0.347 854 845 137 454 0.652 145 154 862 546 ±ξi wi 積分区間の変換 ! _b f (x) dx = ! ₁ f " b − a ξ + b + a # b − a dξ x = b − a₂ ξ + b + a 2 と変数変換すると，

２重積分，３重積分の数値積分

! ₁ −1 ! ₁ −1 ! ₁ −1 f (ξ, η, ζ ) dξdηdζ = n " i=1 n " j=1 n " k=1 f (ξ_i, η_j, ζ_k) w_iw_jw_k ! ₁ −1 ! ₁ −1 f (ξ, η) dξdη = n " i=1 n " j=1 f (ξi, ηj) wiwj

応力の計算

   ε_x ε_y γ_xy    = n ' j=1        ∂N j ∂x 0 0 ∂ N_∂yj ∂N j ∂y ∂N j ∂x        . u j v j / = n ' j=1 0 B j1 . u j v j / １要素の総節点数をnとしたとき，要素内のある点のひずみは次式の用に 形状関数の導関数から計算される． 応力は，これをフックの法則に代入して計算する．すなわち， [D] = E (1 + ν)(1 − 2ν)   1 − ν ν 0 ν 1 − ν 0 0 0 1−2ν₂      σ_x σ_y τ_xy    = [D]    ε_x ε_y γ_xy    = ' D ( n ) j=1 ' B j (* u j v j + ただし，

応力の計算

   σ_x σ_y τ_xy    = [D]    ε_x ε_y γ_xy    = ' D ( n ) j=1 ' B j (* u j v j + 応力を計算する点は，要素積分で用いるガウスポイント（ガウスの数値積 分公式の積分点）を用いる．ガウスポイントにおける形状関数の導関数を 計算して，上式のマトリックス[Bj_{]を計算して各点の応力を計算し，さら} に主応力，応力テンソルの不変量，八面体剪断応力を計算し，要素内で補 間してコンター図などを表示する． ξ η x y x y x y ガウスポイント 実空間での評価点 主応力・主方向 相当応力コンター図

有限要素法のまとめ

•

弱形式，要素分割，変位の有限要素基底による近似，

ガラーキン法の適用

•

弱形式と強形式

•

重み付き残差法と弱形式

•

１次元問題の例

•

２次元問題の例

•

要素座標のパラメータと各種要素

•

ガウスの数値積分公式

•

応力の計算