有限要素法による粘弾塑性体の解析
中央大学 理工学部 教授 川原 睦人1
添字つき記法と加算規約
例えば，2 組の定数 (a1, a2,‥‥‥ an)と (b1, b2,‥‥‥ bn)との積和 A を表わすことを考える． すなわち A = a1b1+ a2b2 +‥‥‥ + anbn (1.1) であるが，(1.1) 式は∑記号を用いて， A = n ∑ i=1 aibi (1.2) と表わすこともできる．よく知られるように n ∑ i=1 という記号は，i という添字にういて，1 から n まで aibiという積をつくり，これを加え合わせることを意味している．しかし，よりいっそう簡 潔にするために n ∑ i=1 記号を省略してしまい， A = aibi (1.3) と書くこともできるであろう．ここで，(1.3) 式のように，同一項に同じ添字がくり返して表われ る場合には，その添字の最大個数まで，それぞれの積をつくりすべてを加え合わせるものと約束 する．すなわち (1.3) 式は (1.1) 式と同じ意味であると解釈する．この約束を加算規約 (Summation convention)と呼んでいる．次に，連立一次方程式を考える． a11x1+ a12x2+‥‥‥ + a1nxn = b1 a21x1+ a22x2+‥‥‥ + a2nxn = b2 (1.4) ‥‥‥‥‥‥ an1x1+ an2x2 +‥‥‥ + annxn= bn (1.4)式は，われわれは，行列によって， [A]· {x} = {b} (1.5) と表わす方法を知っている．しかし，(1.4) 式はまた n ∑ j=1 aijxj = bi(i = 1, 2,· · · n) (1.6)と書くこともできる．ここで，加算規約を用いると (1.6) 式の n ∑ j=1 記号を省略することができ、i は 1から n まで動くことを暗黙のうちに認めて， aijxj = bi (1.7) で表示される．(1.7) 式の意味は (1.4) 式，(1.5) 式の意味とまったく同じものであり，(1.4) 式のよ うに書くと，記述場所を多くとり，長たらしくなるので (1.7) 式のように書くことが行われている． 変位 u, v, w とひずみ εx, εy, εx, γyz, γzx, γxyとの関係は次のように与えられる． εx = ∂u ∂x γyz= 1 2 (∂v ∂z + ∂w ∂y ) εy = ∂v ∂y γzx = 1 2 (∂w ∂x + ∂u ∂z ) (1.8) εz = ∂w ∂z γxy = 1 2 (∂u ∂y + ∂v ∂x ) (1.8)式は添字記法によって， εij = 1 2(ui,j + uj,i) (1.9) と表わすことができる．ここに，i, j は x, y, z をそれぞれ代表しており，, i は∂xi∂ を示すものとする． 応力のつりあい方程式は，物体力を，fx, fy, fzとして， ∂σx ∂x + ∂τxy ∂y + ∂τxz ∂z + fx = 0 ∂τyx ∂x + ∂σy ∂y + ∂τyz ∂z + fy = 0 (1.10) ∂τzx ∂x + ∂τzy ∂y + ∂σz ∂z + fz = 0 (1.10)式で与えられる．これを加算規約によって表示すれば， σij,j + fi = 0 (1.11) ときわめて簡単になる． 弾性体の応力ひずみ方程式は，ラメの定数 λ, µ を用いることによって， σx = (2µ + λ)εx+ λεy+ λεz σy = λεx+ (2µ + λ)εy+ λεz (1.12) σz = λεx+ λεy + (2µ + λ)εz τxy = 2µγxy τyz = 2µγyz τzx= 2µγzx で与えられるが，これはさらに， σij = Eijkl· εkl (1.13) Eijkl= λδijδkl+ µ(δikδjl+ δilδjk) となる．ここに δijはクロネッカーのデルタである． このように加算規約を用いて添字つき記法で表わすと，複雑な式を記述場所を多くとらずに，簡 潔に表示でき，見通しの良い理論展開が可能になる．さらに添字によりその量が明確になってい ることから，電子計算機のプログラミングがきわめて簡単に行うことができる．
2
基礎方程式
有限要素法は，変分原理に基礎をおく，微分方程式の離散化手法である．連続体力学への適用に 際しての基礎方程式を微小変形の仮定のもとに整頓しておく．連続体を任意の大きさをもつ有限 要素と呼ばれる領域に分割し，その体積を V ，表面積を S で表わすことになる． 有限要素 V 内の変位を uiとすれば，ひずみ成分 εij は，ひずみ変位方程式によって，次のよう に与えられる． εij = 1 2(ui,j + uj,i) (2.1) 応力のつりあい方程式は， σij,j + fi = 0 (2.2) である．(2.2) 式の両辺に仮想変位 u∗i をかけ，有限要素全体にわたって積分すると， ∫ V (σij,ju∗i)dV + ∫ V (fiu∗i)dV = 0 (2.3) である．(2.3) 式第 1 項を Gauss の定理によって変形して整頓すれば，結局次の方程式が得られる． ∫ V (σijε∗ij)dV = ∫ S (Piu∗i)dS + ∫ V (fiu∗i)dV (2.4) ここに， ε∗ij = 1 2(u ∗ i,j + u∗j,i) (2.5) Pi = σijnj (2.6) である．njは表面 S にたてた垂線の外向き単位ベクトルの成分である．(2.4) 式が有限要素につい て立てた仮想仕事の方程式である．有限要素の第 α 節点における i 方向の変位を uαiと表わして， 要素内の変位 uiを形状関数 Φαを用いて次のように近似する． ui = Φαuαi (2.7) また，同様に，仮想変位 u∗i についても， u∗i = Φαu∗αi (2.8) が成立するものとする．(2.8) 式を (2.5) 式に代入して整頓すれば， 2ε∗ij = Φα,ju∗αi + Φα,iu∗αj (2.9) であるから，これを (2.4) 式に代入すれば， u∗αi ∫ V (Φα,jσij)dV = u∗αi ∫ S (ΦαPi)dS + u∗αi ∫ V (Φαfi)dV (2.10) (2.10)式を得ることができる． (2.10)式において，u∗αiは任意に選ぶことができることから， ∫ V (Φα,jσij)dV = Ωαi (2.11) Ωαi = ∫ S (ΦαPi)dS + ∫ V (Φαfi)dV (2.12)(2.11)式が成立することが解る．(2.11) 式は，有限要素法の基礎方程式である．Ωαiは，周辺に作 用する表面力 Piおよび要素内部の物体力 fiに形状関数 Φαを重みとしてかけ平均したものになっ ており，換算節点力を表わしている．一方弾塑性解析などでは，応力の増分とひずみの増分との 関係が与えられるので，増分で表わした基礎方程式が必要になる． この場合には，(2.2) 式を増分型式で表わしておき，同様に変形すれば， ∫ V (Φα,j˙σij)dV = ˙Ωαi (2.13) ˙ Ωαi = ∫ S (ΦαP˙i)dS + ∫ V (Φαf˙i)dV (2.14) なる方程式が得られる．
3
弾性体解析
基礎方程式，(2.11) 式あるいは (2.14) 式は，連続体を構成する材料の特性を表わす構成方程式は 用いられずに誘導されている．すなわち，(2.11) 式や (2.14) 式は材料の特性には関係なく成立す る方程式である． 弾性体の構成方程式は，弾性ポテンシャル W によって定義される．すなわち，弾性体とは，次 の方程式を満足する連続体である． σij = ∂W ∂εij (3.1) ここでは，変形は，微小であるとしているから，弾性ポテンシャル W を次のように仮定する． W = 1 2Eijklεijεkl (3.2) ここに Eijklは応力の対称性より 21 個の定数を表わす．等弾性体の場合には，この定数は 2 個に 限定され， Eijkl= λδijδkl+ µ(δikδjl+ δilδjk) (3.3) となる．ここに，λ, µ はラメの定数で，弾性係数，ポアソン比をそれぞれ E, ν とすれば，次のよ うになる． λ = νE (1− 2ν)(1 + ν) µ = E 2(1 + ν) (3.4) (3.2)式を，(3.3) 式に代入すれば， σij = Eijklεkl (3.5) となる．(3.5) 式に (2.7) 式を代入すると σij = Eijkl 1 2(Φβ,luβk + Φβ,kuβl) = EijklΦβ,luβk (3.6) であるから，これを基礎方程式 (2.11) 式に代入して整頓すれば，換算節点力 Ωαiと節点変位 uβkの 関係を次のごとく得ることができる． Kαiβk· uβk = Ωαi (3.7) Kαiβk = ∫ V (Φα,j · Eijkl· Φβl)dV (3.8) Ωαi = ∫ S (ΦαPi)dS + ∫ V (Φαfi)dVKαiβkは一般に剛性行列と呼ばれている． (3.7)式が微小変形の仮定の下における弾性体の解式である．
4
弾塑性体解析
弾塑性体の構成方程式は以下のようにして応力増分 ˙σijとひずみ増分 ˙εij との間の関係で与えら れる． まず，ひずみ成分 εij，弾性ひずみ成分 ε′ij と塑性ひずみ成分 ε′′ijとの和により成り立っているも のと考える. εij = ε′ij + ε′′ij (4.1) 材料の降伏を規定する降伏関数 f を応力 σij，塑性ひずみ ε′′ij の関数として与えることにし，ひず み硬化係数 π を用いて， f (σij, ε′′ij) = π (4.2) と表わす．降伏関数を用いて，それぞれ次の状態を定義する． i）初期弾性 ( ˙ε′′ij = 0) f < π, ˙π = 0 (4.3) ii）負荷 (ε′′ij ̸= 0) f = π, ˙π̸= 0, ∂f ∂σij ˙σij > 0 (4.4) ii）除荷および中立負荷 ( ˙ε′′ij = 0) f = π, ˙π = 0, ∂f ∂σij σij ≤ 0 (4.5) (4.5)式第 3 式の等号が成立するとき中立負荷と呼ばれている．(4.4) 式が成立するときのみ塑性 ひずみ成分の増分 ˙ε′′ij が生じ，これが流れ法則に従うものと考える．すなわち， ˙ ε′′ij = Λ ∂f ∂σij (4.6) 弾性ひずみ増分については，(3.5) 式が成立するとし，増分で表すと， ˙σij = Eijkl· ˙ε′kl (4.7) となる．また，(4.2) 式における π は ˙ε′′ij のみの関数であるとして， ∂f ∂σij · ˙σ ij+ ∂f ∂ε′′ij · ˙ε ′′ ij = ∂π ∂ε′′ij · ˙ε ′′ ij (4.8) なる関係が得られる．(4.1)，(4.6)，(4.7)，(4.8) 式より Λ を求め，これをあらためて (4.1)，(4.6)， (4.7)式に代入して整頓すると結局次のように整頓することができる． ˙σij = Cijkl· ˙εij (4.9) Cijkl = Eijkl− Eijpq ( ∂f ∂σpq ) ·( ∂f ∂σrs ) Erskl ( ∂f ∂σmn )( ∂π ∂ε′′ij − ∂f ∂ε′′mn + Epqmn ∂f ∂σpq )(4.9)式は (4.4) 式が成立するときに用いられる． よって弾塑性体の構成方程式を整頓すると次のようになる． ˙σij = Dijkl· ˙εkl i）初期弾性 f < π, ˙π = 0のとき Dijkl = Eijkl ii）負荷 f = π, ˙π ̸= 0, ∂f ∂σij ˙σij > 0のとき Dijkl = Cijkl iii）除荷および中立負荷 f = π, ˙π = 0, ∂f ∂σij ˙σij ≤ 0 のとき Dijkl = Cijkl (4.10)式に (2.7) 式を代入すれば， ˙σij = DijklΦβ,l˙uβk (4.10) となるから，(2.13) 式に代入して換算節点力の増分と，節点変位の増分との関係を次のごとく得 ることができる． KαiβkUβk = Ω˙αi (4.11) Kαiβk = ∫ V (Φα,j· Dijkl· Φβ,l)dV ˙ Ωαi = ∫ V S (ΦαP˙i)dS + ∫ V (Φαfi)dV Dijklはその要素が塑性状態にあるとき，基準となる σijと ε′′ijなどの関数で増分 ˙σij, ˙εijなどは含ん でいない．
5
粘弾性体解析
粘弾性体の構成方程式は緩和関数などによる積分表示される場合とレオロジーモデルを用いた 微分表示で与えられる場合がある． 緩和関数を用いる場合には，構成方程式は次のようになる． σij = ∫ t 0 Gijkl(t− τ) · ∂εkl ∂τ · dτ (5.1) ここに，Gijklは緩和関数で，多要素 M axwell 型の場合には， Gijkl = G0δijδkl+ (δikδjl+ δilδjk) N ∑ n=1 Gnexp (t− τ τn ) (5.2) が用いられる．Gnは n 番目の要素の緩和関数，τnは緩和時間である． (2.7)式より ∂εij ∂τ = 1 2 ( Φα,j ∂Uαi ∂τ + Φα,i ∂Uαj ∂τ ) (5.3)であるから，これを (5.1) 式に代入して整頓すれば，次の関係を得ることができる． σij = ∫ t 0 Gijkl(t− τ) ∂Uβl ∂τ dτ· Φβ,k (5.4) いま，時間を適当に分割し，その分割された時間区間ごとに解を追跡する逐次近似法によって解 析することにする． 着目する分割された時間区間をあらためて 0∼t 区間とし，この区間内の変位を次のごとく近似 的に二次式で与えられるとする． Ui(τ ) = t2− τ2 t2 Ui(0) + τ2 t2Ui(t) + τ (t− τ) t ˙ Ui(0) (5.5) これを時間で微分すれば， ∂Ui(t) ∂τ = t2− 2τ t2 Ui(0) + 2τ t2Ui(t) + t− 2τ t ˙ Ui(0) (5.6) となる．これらは節点変位についても成立するから (5.4) 式に代入して整頓すれば，結局次のよう に整頓することができる． σij = Hijkl· Φβ,lUβk− σij∗ (5.7) Hijkl = ∫ t 0 2Gijkl(t− τ) · τ2 t2 · dτ σij∗ = ∫ t 0 Gijkl(t− τ) t2− 2τ t2 dτ· Φβ,kUβl(0) + ∫ t 0 Gijkl(t− τ) t− 2τ t dτ · Φβ,kUβl(0) (5.9)式を基礎方程式 (2.11) 式に代入すると，
Kαiβk· Uβk = Ωαi + Ω∗αi (5.8)
Kαiβk = ∫ V (Φα,j · Hijkl· Φβ,l)dV Ω∗αi = ∫ V (Φα,jσ∗ij)dV のごとく，0∼t 区間に対する解式を作成することができる． 粘弾性体の構成方程式は，(5.1) 式で与えられると同時に，微分表示することも可能である． すなわち， Qαij = Aαβijkl· qβkl+ Bijklαβ + ˙gklβ (5.9) i, j, k, l = 1, 2, 3 α, β = 1, 2,· · · L ここに， Q1ij = σij, qkl1 = εkl Qαij = 0(α̸= 1), qklβ = hβkl(β ̸= 1) である．L はレオロジーモデルの要素の総数に対応し，hβkl(β ̸= 1) は潜在変数と呼ばれ，直接測定 することはできない量である．レオロジーモデルの連結の仕方を表わす定数の組 Aαβijkl, B αβ ijklは機
械的に作成することができる．さて，解を求めようとする時間を，適当な時間間隔に分割し，その うちの一つを ∆t と表わすことにする．∆t 区間において，ひずみ速度 ˙εij が一定である場合には， ˙ qαij = q α ij − qijα(0) ∆t (5.10) と近似することができる．(5.10) 式と (5.9) 式に代入して整頓し，ひずみと潜在変数とを分けて表 わせば， σij = K (1) ijkl· εkl+ K (2)α ijkl · h α kl + Hijkl(1) · εkl(0) + H (2)α ijkl · h α kl(0) (5.11) Kijkl(2)α· εij + K (3)αβ ijkl · h β ij + Hijkl(2)α· εij(0) + H (3)αβ ijkl · h β ij(0) = 0 (5.12) である． Kijkl(1) = A(1)ijkl+ 1 ∆tB (1) ijkl, Hijkl =− 1 ∆tB (1) ijkl, Kijkl(2)α = A(2)αijkl + 1 ∆tB (2)α ijkl, H α ijkl =− 1 ∆tB (2)α ijkl, Kijkl(3)αβ = A(3)ijkl+ 1 ∆tB (3)αβ ijkl , H αβ ijkl =− 1 ∆tB (3)αβ ijkl , である．(5.12) 式を用いて (5.11) 式から hβijを消去すれば，結局次のように構成方程式を整頓す ることができる． σij = Dijkl· εkl− σij∗ (5.13) hβij = −Sβ(1)· εij − S (2) β · εij(0)− S (3) αβ · h β ij(0)
Sβ(1) = Kijkl−1(3)αβ· Kijkl(2)αβ, Sβ(2) = Kijkl−1(3)αβ· Hijkl(2)α, Sαβ(3) = Kijkl−1(3)αγ · Hijkl(3)γα, Dijkl = Kijkl(1) + K (2)γ ijkl · K −1(3)γδ ijkl · K (2)δ ijkl, σij∗ = (Hijkl(1) − Kijkl(2)γ· Sγ(3))· εkl(0) + ( Hijkl(2)γ− Kijkl(2)δ· Sδγ(3) ) · hγ kl(0), (5.13)式における Dijklは弾性体の応力ひずみ関係を表わす係数と相似な関係にあり，粘弾性定数 と時間区間 ∆t とにより成り立っている．また σ∗ijは ∆t 区間の初期値により構成されており，既 知項である．(5.12) 式を用いれば，ひずみ εij より潜在変数 hαijを計算することができる． (5.13)式を (2.11) 式に代入すれば，粘弾性体の解式を次のように得ることができる．
Kαiβk· Uβk = Ωαi+ Ω∗αi (5.14)
ここに， Kαi = ∫ V (Φα,j · Dijkl· Φβ,k)dV Ω∗αi = ∫ V (Ωα,jσ∗ij)dV である．(5.14) 式は (5.8) 式と形式的にはまったく同一の方程式となる．さらに，この方程式は， 初期応力をもつ弾性体の方程式と相似になり，粘弾性体の解析は，時間間隔ごとに弾性体と同じ 計算と逐次行うことによって進められる．
6
粘弾塑性体解析
粘弾塑性体とは，ある降伏条件が満足されないうちは弾性体として挙動し，降伏条件が満足さ れた後には，粘塑性的に変形する物体のことを言う．この特質をもつ材料の構成方程式としては， Praqerが提案する方程式が代表的である．等方性の材料を対象とすることにして，まず偏差応力 Sijと偏差ひずみ eij を導入する． Sij = σij− σkk 3 δij (6.1) eij = εij − εkk 3 δij (6.2) 降伏条件は，弾塑性体の場合と同様にして f (σij, e′′ij = π (6.3) で与えられるものとする．ここに，e′′ij は粘弾性偏差ひずみを表わしている．また，偏差ひずみ eij は弾性偏差ひずみ e′ijと粘弾性偏差ひずみ e′′ij の和として与えられるとする． eij = e′ij + e′′ij (6.4) 一方体積ひずみ ε と等圧応力 σ はそれぞれ次のように定義されるが， σ = σkk 3 (6.5) ε = εkk 3 (6.6) これ等は，弾性成分のみより成立すると仮定する．すなわち，粘弾塑性ひずみは偏差ひずみにし か起こらず，体積変形は，すべての弾性変形としてのみ挙動するものとする． 以上のことから， σ = 3K · ε (6.7) Sij = 2G· e′ij (6.8) である．ここに，K は体積弾性係数を G はせん断弾性係数をそれぞれ表わしている．(6.8) 式を変 形して，時間に対する増分の型式で表わせば， ˙Sij = 2G( ˙eij − ˙e′′ij) = 2G( ˙eij − ˙ε′′ij) (6.9)
となる．ここに，ε′′ijはひずみの粘弾塑性成分で，弾性成分 ε′ijにより， εij = ε′ij + ε′′ij (6.10) と表わすことができる．さらに， ε′′ii= 0 (6.11) であるから，(6.9) 式が成立する． 一方 Praqer によれば，粘弾塑性ひずみ速度 ˙ε′′ij は次のように偏差応力で表わすことができると されている． ˙ ε′′ij = 1 2η(1− kJ −1 2 2 )· Sij (6.12)
ここに，J2は J2 = 1 2SijSij (6.13) のごとく偏差応力の 2 次式の不変量であり，η と k は実験により定まる定数である． (6.7)式と (6.9) 式により ˙σij ={µ(δikδjl+ δliδkj) + λδijδkl} ˙εkl − µ(δikδlj + δliδkj) ˙ε′′kl (6.14) λ = 1 3(3K− 2G) = νE (1 + ν)(1− 2ν) µ = G = E 2(1 + ν) となる．(6.14), (6.12), (6.13) 式によって，(6.3) 式が成り立つ場合の構成方程式が与えられ，この とき (6.11) 式が付帯条件となる． (6.13)式が満足されないときは，通常の弾性体と同じ構成方程式によって表わすことができる． 以上をまとめて次のごとく表わす． ˙σij = Eijkl· ˙εkl− σ∗ij (6.15) σij∗ = 2µ· ˙ε′′ij ˙ ε′′ij = 1 2η(1− kJ −1 2 2 )· Sij, J2 = 1 2SijSij Eijkl = µ(δikδjl+ δilδjk) + λδijδkl (2.7)式を (6.15) 式に代入し，その方程式を (2.11) 式に代入して整頓すれば，最終的に有限要素法 の解式として，次の方程式を得ることができる．
Kαiβk· ˙Uβk = Ω˙αi+ Ω∗αi (6.16)
kαiβk = ∫ V (Φα,j· Eijkl· Φβ,l)dV Ω∗αi = ∫ V (Φα,j· σij∗)dV (6.16)式は，時間に関する一階の微分方程式であるので，これを時間間隔を細かく分割して逐次 積分を進める方法を用いることになる． (6.16)式左辺係数 Kαiβkはいまの場合弾性体の係数と同一のものとなり， ˙Ω∗αiを修正することで 計算が進行する．
7
熱連成弾性体解析
この節では，熱の影響を受ける弾性体の解析を扱うことにする．まず，基礎方程式はひずみ変位 方程式 (2.1) 式，応力つりあい方程式 (2.2) 式，の外に温度に関する方程式と，構成方程式とが必 要になる． 温度勾配 ζiは相対温度 θ の勾配として与えることができ， ζi = θ,i (7.1)となる．エネルギーの平衝方程式は，ε を単位質量当りの内部エネルギーとして， ρ ˙ε = σijε˙ij + qi,i+ ρh (7.2) となる．ここに ρ は密度，qiは熱流束，h は単位質量当りの熱発生である．単位質量当りのエント ロピーを η で表わすと，Clasius-Duhem の不等式は， ρR ˙η− qi,i− ρh + 1 TqiT,i ≥ 0 (7.3) となる．ここに T は絶対温度である． ここで，自由エネルギー φ，および内部減衰 σ をそれぞれ次のように導入する． φ = ε− ηT (7.4) σ = σijε− ρ( ˙φ + η ˙T ) (7.5) (7.4), (7.5)式を用いて (7.2) 式を変形すると次のごとくになる． ρ ˙φ = σijεij − ρη ˙T − σ (7.6) ρT ˙η = qi,i+ ρh + σ (7.7) (7.6)式は構成方程式の誘導に用いられ，(7.7) 式はエントロピーの平衝方程式である． (7.7)式の両辺に θ ∗ を掛け，線型化するために右辺第 3 項を省略して，変形し表面 S，体積 V をもつ領域について積分し整頓すると，結局次の方程式が得られる． ∫ V (θ∗ρT0η)dV +˙ ∫ V (θ,i∗qi)dV = ∫ S (θ∗qiηi)dS + ∫ V (θ∗ρh)dV (7.8) ここに，絶対温度 T は，自然状態の温度 T0に近似的に等しいとする． (2.7)式と同様にして，温度 θ も，要素の各節点における温度 θαを用いて，形状関数 Φαによっ て，次のように記述されるとする． θ = Φαθα (7.9) (7.9)式を (7.8) 式に代入して整頓し，α 節点での値 θα∗ は任意であることを用いると，(2.11) 式に 対応する方程式として， ∫ V (ΦαρT0η)dV +˙ ∫ V (Φαqini)dV = Γα (7.10) Γα = ∫ S (Φαqini)dS + ∫ V (Φαρh)dV (7.11) なる温度に関する基礎方程式が得られる． 構成方程式は，自由エネルギー φ が次のように与えられるとして構成される． φ = 1 2Eijklεijεkl− Bijεijθ− ρC 2T0 (7.12) これより， σij = ∂φ ∂εij = Eijklεkl− Bijθ (7.13) ρη = −∂φ ∂θ = Bijεij + ρC T0 θ (7.14)
である．ここに， Eijkl = λδijδkl+ µ(δikδjl+ δilδjk) Bij = b· δij b = E 1− 2να となる． λ，µ はラメの定数，α は線膨張係数，C は比熱である．一方熱流束 qi は，フーリエの法則に よって， qi = κζi (7.15) と与えられる．κ は熱伝導率である． (7.13), (7.14), (7.15)式にそれぞれ，(2.7), (7.9) 式を代入して整頓すると，次のように方程式を 整頓することができる． σij = EijklΦβ,lUβk − BijΦβ· θβ (7.16) ρη = BijΦβ,iUβj + ρC T0 Φβθβ (7.17) qi = κΦβ,iθβ (7.18) (7.16)式を (2.11) 式に代入して整頓すると，応力変形に関する方程式が得られる．
Kαiβk· Uβk− Tαiβ· θ = Ωαi (7.19)
ここに， Kαiβk = ∫ V (Φα,jEijklΦβ,l)dV Tαiβ = ∫ V (Φα,jBijΦβ)dV Ωαi = ∫ S (Φαρi)dS + ∫ V (Φαfi)dV である． 一方，(7.17) 式，(7.18) 式を (7.10) 式に代入して整頓すれば，結局，熱伝導に関する方程式を (7.20)式のように得ることができる． Mαβθ˙β + NαiβU˙αi+ Jαβθβ = Γα (7.20) Mαβ = ∫ V (ΦαρCΦβ)dV Nαiβ = ∫ V (ΦαT BijΦβ,j)dV Jαβ = ∫ V (Φα,i· κ · Φβ,i)dV Γα = ∫ V (Φαρh)dV + ∫ S (Φαqini)dS (7.19)式と (7.20) 式を連立して解けば，熱連成弾性体の解析を行うことができる．(7.20) 式には， 変位 Uαiが連成しているが，この項を省略して， Jαβθ = Γα (7.21) とすれば，通常の熱伝導の方程式を得ることになる．
8
有限変形弾性体解析
直角座標系 xiによる物質表現を用いいることにする．変形前後の物質座標をそれぞれ xi, ziと表 わせば変位 Uiは次のように与えられることになる． Ui = zi− xi (8.1) これにより変形勾配 Fijは， Fijzi,j = δij + Ui,j (8.2) となる．Green のひずみ εij は変位 Uiによって (8.3) 式のごとく表わすことができる． εij = 1 2(Ui,j + Uj,i+ Um,jUm,j) (8.3) Kirchhoffの応力を Sijとすれば，応力の平衝方程式は，ρ を密度，fiを物体力として， (FijSjk),k+ ρfi = 0 (8.4) で与えられる．ここに慣性力は無視されている．仮想仕事の方程式を誘導するために (8.4) 式の両 辺に Ui∗をかけて，変形前の体積 V0，表面 S0をもつ領域全体について積分すると次の方程式が得 られる． ∫ V0 { (FijSjk),kUi∗ } dV ∫ V0 (ρ0fiUi∗)dV = 0 (8.5) (8.5)式第 1 項に Green-Gauss の定理を用いて変形し，整頓すると，結局次の方程式を得ることが できる． ∫ V0 (Sijε∗ij)dV = Ω (8.6) ε∗ij = 1 2(U ∗i,j + Uj,i∗ + Um,i∗ Um,j+ Um,iUm,j∗ ) (8.7)
Ω = ∫ S0 [ Sijnj(δki+ Uk,i)Uk∗ ] dS + ∫ V0 (ρ0fiUi∗)dV ここに njは表面 A0にたてた単位法線ベクトルの成分である．構成方程式としては，簡単のため に，線型弾性体と同じものを用いることにする． Sij = Eijkl· εkl (8.8) さて，解析する物体を有限の大きさをもつ有限要素に分割し，その有限要素の第 α 節点第 i 方向 変位を Uαiと表わす．要素内の変位 Uiを変位の形状関数 Φαによって， Ui = ΦαUαi (8.9) と近似する．(8.9) 式を (8.3) 式，(8.7) 式に代入すれば，
2· εij = Φα,iUαj + Φα,jUαi + Φα,iΦβ,jUαmUβm (8.10)
2εij = Φα,iUαj∗ + Φα,jUαi∗ + Φα,iΦβ,jUαm∗ Uβm+ Φα,iΦβ,jUαmUβm∗ (8.11)
である．(8.8) 式を (8.6) 式に代入して整頓すれば，
∫
V0
であるから，これに，(8.10) 式，(8.11) 式を代入して整頓すれば，有限変形弾性体の解式を次のよ うに得ることができる． KαiβjUβj = Ωαi (8.13) Kαiβj = ∫ V0 [ Φα,k(δli+ Φδ,lUδi)Eklmn· (δmj + 1 2Φγ,mUγj)Φβ,n ] dV0 Ωαi = ∫ S0 (ρiΦα)dS + ∫ V0 (ρ0fiΦα)dV Pk = Sij(δki+ Φβ,iUβk)nj (8.13)式を Newton-Raphson 法などにより解けば，解を得ることができる．
