主不変量を用いた Ogden 材料モデルの定式化

4.2.1 主不変量を用いたOgden材料モデル^(4.7)

Ogden材料モデルは，非圧縮性ゴム状物質に対する等方性弾性ポテンシャルW として，

次式のように定義される．

W = XN

r=1

µ_r

α_r(λ^α₁^r+λ^α₂^r +λ^α₃^r −3) (4.1)

λ₁λ₂λ₃= 1 (4.2)

ただし，N は正の定数で，µ_rおよびα_rは，µ_rα_r>0(k= 1,2,· · · , N)を満足する材料定 数である．ここで，λ_iは右ストレッチテンソルU(=C^1/2，ただしCは右Cauchy-Green 変形テンソル)の主値である．上式で，N = 2，µ₁ = 2c₁,α₁ = 2,µ₂ =−2c₂,α₂ =−2の とき，次式のMooney-Rivlin材料モデルの式に帰着する．

W =c₁(I−3) +c₂(II−3) (4.3)

III = 1 (4.4)

ただし，c₁，c₂はMooney定数，I，II，IIIはCの第1，2，3主不変量である．式(4.2)，

(4.4)は非圧縮条件を示す．ここで，後の式展開の利便性を考え，式(4.1)，(4.2)を次式の

ようにCの主値Λ_i=λ²_i を用いて書き換える．

W = XN

r=1

µ_r

α_r(Λ^α₁^r/2+ Λ^α₂^r/2+ Λ^α₃^r/2−3) (4.5)

Λ₁Λ₂Λ₃ = 1 (4.6)

上のOgden式は，Cの主値を用いるため，超弾性体の定義にしたがって応力テンソル，構

成則テンソルを求める際には，固有値解析を全ての積分点において実施し，さらに基底ベ クトルを主軸方向へ変換させる数値的な操作を行わなければならない．

Basar^(4.7)らは，次式のように，Λ_iを主不変量I，II，IIIの関数として表した．

Λ_i = Λ_i(I, II, III) (4.7)

上式から，式(4.5)のOgden材料モデルの構成則を主不変量の関数として書き直せる．一

般に，Ogden材料モデルの応力テンソル・構成則テンソルは，主値に対応した固有ベクト

ルを基底とした座標系に変換してから求める必要がある．以上のように構成則をあらかじ め主不変量を用いて表しておけば，Mooney-Rivlin材料モデルのときと同じように主不変 量の微分計算を考えればよく，Cから主値を求める固有値解析と固有ベクトルを基底した 座標系への座標変換という数値的な操作を省くことができる．

まず，超弾性体の定義より，第2種Piola-Kirchhoff応力テンソルS，構成則テンソルD はそれぞれ次式のように与えられる．

S = 2∂W

∂C (4.8)

D= 4 ∂²W

∂C∂C (4.9)

W はΛ_iの関数であるから，微分の連鎖則より，式(4.8)，(4.9)はそれぞれ次式のように なる．

S = 2∂W

∂Λ_i

∂Λ_i

∂C (4.10)

D = 4 ∂²W

∂Λ_i∂Λ_j

∂Λ_i

∂C ⊗∂Λ_j

∂C + 4∂W

∂Λ_i

∂²Λ_i

∂C∂C (4.11)

ここで，式(4.7)より，Wは主不変量I，II，IIIの関数でもあることから，式(4.10)，(4.11) を次式のように書き直すことを考える．

S = 2 µ

a_I ∂I

∂C +a_II∂II

∂C +a_III∂III

∂C

¶

(4.12)

D = 4 µ

b_I ∂²I

∂C∂C +b_II ∂²II

∂C∂C +b_III ∂²III

∂C∂C +c_{I I} ∂I

∂C ⊗ ∂I

∂C +c_{II II}∂II

∂C ⊗∂II

∂C +c_{III III}∂III

∂C ⊗ ∂III

∂C +c_{I II} ∂I

∂C ⊗∂II

∂C +c_{I III} ∂I

∂C ⊗∂III

∂C +c_{II III}∂II

∂C ⊗∂III

∂C

¶

(4.13) ただし，a_K，b_K，c_{K L}(K, L=I, II, III)は式(4.10)，(4.11)から式(4.12)，(4.13)へ変換す る際に求められる未知係数である．ここで，a_I =b_I =c₁，a_II =b_II =c₂，a_III =b_III = 0， c_{K L} = 0のとき，式(4.12)，(4.13)はMooney-Rivlin材料モデルに帰着する．

4.2.2 未知係数a_K，b_K，c_{K L}の求め方^(4.9)

式(4.7)は，次式の特性方程式の根として導かれる．

Λ³−IΛ²+IIΛ−III = 0 (4.14)

Λ_iの条件Λ₁ 6= Λ₂ 6= Λ₃，Λ₁ >Λ₂ = Λ₃，Λ₁ <Λ₂ = Λ₃，Λ₁ = Λ₂ = Λ₃ によって，式

(4.14)の解が異なるため，それぞれの場合について以下に示す．

いま，次の4つの関数を定義する．

f(I, II, III) =I²−3II (4.15)

g(I, II, III) =27III²+ 4II³I²II²+ 4I³III−18IIIIII (4.16) h(I, II, III) =2I³

27 −III

3 +III (4.17)

w(I, II, III) =1 3cos⁻¹

·2I³−9III+ 27III 2(I²−3II)^3/2

¸

(4.18) このとき，上の4つのΛ_iの条件について，以下のように場合分けを行う．

(1) f = 0のとき，Λ₁ = Λ₂= Λ₃

Λ₁ = Λ₂ = Λ₃= 1

3I (4.19)

(2) f 6= 0，g= 0，h >0のとき，Λ₁>Λ₂ = Λ₃ Λ₁= 1

3 h

I+ 2p

I²−3II i

(4.20) Λ₂= Λ₃ = 1

3 h

I−p

I²−3II i

(4.21) (3) f 6= 0，g= 0，h <0のとき，Λ₁<Λ₂ = Λ₃

Λ₁= 1 3

h

I−2p

I²−3II i

(4.22) Λ₂= Λ₃ = 1

3 h

I+p

I²−3II i

(4.23) (4) f 6= 0，g6= 0のとき，Λ₁6= Λ₂6= Λ₃

Λ₁ = 1 3I+2

3

pI²−3IIcosw (4.24)

Λ₂ = 1 3I+2

3

pI²−3IIcos(2

3π−w) (4.25)

Λ₃ = 1 3I+2

3

pI²−3IIcos(2

3π+w) (4.26)

式(4.19)〜(4.26)より，上の4つの条件それぞれに対して，(∂Λ_i)/(∂C)，(∂²Λ_i)/(∂C∂C) と(∂K)/(∂C)，(∂²K)/(∂C∂C)，(K =I, II, III)の間の関係を導き，式(4.10)，(4.11) へ代入すると，未知係数a_K，b_K，c_{K L}を決定することができる．それぞれの場合のa_K， b_K，c_{K L}を以下に示す．なお，具体的な計算方法は文献^(4.7)の通りである．

まず，次の関数を定義する．

¯

ω(Λ) =X

r

µ_r

α_r(Λ^αr/2−1), W = ¯ω(Λ₁) + ¯ω(Λ₂) + ¯ω(Λ₃)

¯

ω⁰(Λ) = d¯ω dΛ =X

r

µ_r

2 Λ^(αr/2)−1,

¯

ω⁰⁰(Λ) = d²ω¯ dΛ² =X

r

µ_r 2

³α_r 2 −1

´

Λ^(αr/2)−2 (4.27)

D₁ = (Λ₁−Λ₂)(Λ₁−Λ₃), D₂ = (Λ₂−Λ₃)(Λ₂−Λ₁),

D₃ = (Λ₃−Λ₁)(Λ₃−Λ₂) (4.28)

Ω₁ = ¯ω⁰(Λ₁)−ω¯⁰(Λ₂)−(Λ₁−Λ₂)¯ω⁰⁰(Λ₂),

Ω₂ = ¯ω⁰⁰(Λ₁)−ω¯⁰⁰(Λ₂) (4.29) このとき，a_K，b_K，c_{K L}はΛ_iの条件により，以下のように表される．

(1) Λ₁ = Λ₂ = Λ₃のとき

a_I= ¯ω⁰(Λ), a_II =a_III = 0,

b_I = ¯ω⁰(Λ) + 2Λ¯ω⁰⁰(Λ), b_II =−¯ω⁰⁰(Λ), b_III = 0, c_{I I} = ¯ω⁰⁰(Λ),

c_{II II} =c_{III III} =c_{I II} =c_{I III} =c_{II III} = 0 (4.30)

(2) Λ₁ >Λ₂ = Λ₃またはΛ₁<Λ₂= Λ₃のとき a_I = [¯ω⁰(Λ₁)−ω¯⁰(Λ₂)] Λ²₁

(Λ₁−Λ₂)² + ¯ω⁰(Λ₂) a_II =−[¯ω⁰(Λ₁)−ω¯⁰(Λ₂)] Λ₁

(Λ₁−Λ₂)² a_III = [¯ω⁰(Λ₁)−ω¯⁰(Λ₂)] Λ₁

(Λ₁−Λ₂)²

b_I = Ω₁ Λ²₁

(Λ₁−Λ₂)² + ¯ω⁰(Λ₂) + ¯ω⁰⁰(Λ₂)(Λ₁+ Λ₂) b_II =−Ω₁ Λ₁

(Λ₁−Λ₂)² −ω¯⁰⁰(Λ₂) b_III = Ω₁ 1

(Λ₁−Λ₂)² c_{I I} = Ω₁ −4Λ³₁Λ₂

(Λ₁−Λ₂)⁵ + Ω₂ Λ⁴₁

(Λ₁−Λ₂)⁴ + ¯ω⁰⁰(Λ₂) c_{II II} = Ω₁−2Λ₁(Λ₁+ Λ₂)

(Λ₁−Λ₂)⁵ + Ω₂ Λ²₁ (Λ₁−Λ₂)⁴ c_{III III} = Ω₁ −4

(Λ₁−Λ₂)⁵ + Ω₂ 1 (Λ₁−Λ₂)⁴ c_{I II} = Ω₁Λ²₁(Λ₁+ 3Λ₂)

(Λ₁−Λ₂)⁵ + Ω₂ −Λ³₁ (Λ₁−Λ₂)⁴ c_{I III} =c_{II II}

c_{II III} = Ω₁(3Λ₁+ Λ₂)

(Λ₁−Λ₂)⁵ + Ω₂ −Λ³₁

(Λ₁−Λ₂)⁴ (4.31)

(3) Λ₁ 6= Λ₂ 6= Λ₃のとき a_I =b_I =

X3

i=1

¯

ω⁰(Λ_i)Λ²_i D_i , a_II =b_II =−

X3

i=1

¯ ω⁰(Λ_i)Λ_i

D_i , a_III =b_III =

X3

i=1

¯ ω⁰(Λ_i)

D_i c_{I I} =

X3

i=1

Λ²_i D_i²

·

¯

ω⁰⁰(Λ_i)Λ²_i + 2ω¯⁰(Λ_i)(3III−Λ_iII) D_i

¸

c_{II II} = X3

i=1

1 D²_i

·

¯

ω⁰⁰(Λ_i)Λ²_i + 2ω¯⁰(Λ_i)(III−Λ³_i) D_i

¸

c_{III III} = X3

i=1

1 D²_i

·

¯

ω⁰⁰(Λ_i) + 2ω¯⁰(Λ_i)(I −3Λ_i) D_i

¸

c_{I II} = X3

i=1

Λ_i D²_i

·

−¯ω⁰⁰(Λ_i)Λ²_i +ω¯⁰(Λ_i)(Λ²_iI−3III) D_i

¸

c_{I III} =c_{II II} c_{II III} =

X3

i=1

1 D_i²

·

−¯ω⁰⁰(Λ_i)Λ_i+ ω¯⁰(Λ_i)(3Λ²_iI−II) D_i

¸

(4.32)

なお，一般にa_K，b_Kは，Λ_iが2重根，3重根のときはa_K 6=b_Kとなる．また,上の(1)

〜(4)の場合分けの計算を実際に行う際は，式f = 0，g= 0の判定基準として，それぞれ

|f|<1.0×10⁻¹²，|g|<1.0×10⁻⁶とおいたとき，最もゼロ割，オーバーフローなどの計 算実行エラーを回避する傾向にあった．

4.2.3 混合型ソリッド要素の定式化

以上の主不変量を用いたOgden材料モデルの構成則を等容変形，体積変形に分離した 形で混合型ソリッド要素へ適用した例は筆者らの知る限りない．そこで本節では，主不変 量を用いたOgden材料モデルの混合型ソリッド要素への実装方法を示す．

混合法の場合，新たな未知数として不定静水圧pを導入し，次式のように等容変形，

体積変形に分離して第2 種Piola-Kirchhoff 応力テンソル，構成則テンソルを定義する

(4.11)(4.12)(4.13)．

S = 2∂W¯

∂C¯ : ∂C¯

∂C + 2p∂III

∂C (4.33)

D= 4 ¯D+ 4p ∂²III

∂C∂C (4.34)

ここで，W¯，C¯は，それぞれひずみエネルギーW，右Cauchy-Green変形テンソルCか ら体積変形を除いたW，Cの偏差部分である．このとき，W¯ はC¯の関数となり，C¯は次 式のように定義される．

C¯ =III⁻¹³C (4.35)

また，式(4.34)中のD¯ は(∂²W¯)/(∂C¯∂C)¯ のテンソル値テンソル関数で，次式のように 与えられる．

D¯ = ¯D_DEV −2

3S¯⊗C⁻¹−2

3C⁻¹⊗S¯ +2

3J⁻²³ µ

2∂W¯

∂C¯ :C

¶ ·

I_C−1 −1

3C⁻¹⊗C⁻¹

¸

(4.36) ただし，以下に定義されるテンソルを用いる．

S¯ = 2∂W¯

∂C¯ : ∂C¯

∂C (4.37)

(I_C−1)_ijkl ≡ 1 2 h

C_ik⁻¹C_jl⁻¹+C_il⁻¹C_jk⁻¹ i

(4.38) D¯_DEV = 4J⁻⁴³

· ∂²W¯

∂C∂¯ C¯ +1

9(C : ∂²W¯

∂C¯∂C¯ :C)C⁻¹⊗C⁻¹

−1

3C⁻¹⊗( ∂²W¯

∂C¯∂C¯ :C)−1 3( ∂²W¯

∂C¯∂C¯ :C)⊗C⁻¹

¸

(4.39)

以上より，式(4.5)のOgden体の第2種Piola-Kirchhoff応力テンソル，構成則テンソル を主不変量を用いて導く場合は，それぞれ次式の計算が必要となる．

∂W¯

∂C¯ = ¯a_I ∂I¯

∂C¯ + ¯a_II∂II¯

∂C¯ + ¯a_III∂III¯

∂C¯ (4.40)

∂²W¯

∂C∂¯ C¯ = ¯b_I ∂²I¯

∂C∂¯ C¯ + ¯b_II ∂²II¯

∂C¯∂C¯ + ¯b_III ∂²III¯

∂C∂¯ C¯ + ¯c_{I I} ∂I¯

∂C¯ ⊗ ∂I¯

∂C¯ + ¯c_{II II}∂II¯

∂C¯ ⊗∂II¯

∂C¯ + ¯c_{III III}∂III¯

∂C¯ ⊗ ∂III¯

∂C¯ + ¯c_{I II} ∂I¯

∂C¯ ⊗ ∂II¯

∂C¯ + ¯c_{I III} ∂I¯

∂C¯ ⊗∂III¯

∂C¯ + ¯c_{II III}∂II¯

∂C¯ ⊗∂III¯

∂C¯ 

(4.41) ここで，I¯，II¯，III¯ はC¯の第1，2，3主不変量である．また，¯a_K，¯b_K，c¯_{K L}は，前節で 導いた未知係数a_K，b_K，c_{K L}の計算で，Λ_iの替わりにC¯ の主値Λ¯_iを用いることで求め られる．なお，式(4.33)〜(4.41)の中で，以下のテンソルの計算を用いる．

∂I

∂C =I, ∂II

∂C =II −C, ∂III

∂C =IIIC⁻¹,

∂I¯

∂C¯ =I, ∂II¯

∂C¯ = ¯II −C,¯ ∂III¯

∂C¯ = ¯IIIC¯⁻¹,

∂C¯

∂C =III⁻¹³

· I−1

3C⊗C⁻¹

¸

(4.42)

4.2.4 非圧縮超弾性シェル要素の定式化

ここでは，主不変量を用いたOgden材料モデルの構成則を，著者らが開発した新たな 非圧縮超弾性シェル要素^(4.15)へ適用した例を示す．文献^(4.15)のシェル要素は，既往の大 ひずみシェル要素^(4.7)の定式化とは異なり，板厚方向にひずみ仮定場を設けることで，面 外せん断ひずみを考慮しながら板厚を陰的かつ効率的に取り扱っている．本シェル要素で は，構成則の導出に平面応力条件と非圧縮条件を用いる．非圧縮性から生じる不定静水圧 pは，シェルの平面応力条件により消去される．

まず，非圧縮条件III = 1より，Cの板厚方向の成分C₃₃を消去する．主不変量I，II， IIIは，成分C_ij を用いると

I =C_ii (4.43)

II= 1

2(C_iiC_jj −C_ijC_ji) (4.44) III =ε_ijkC_1iC_2jC_3k (4.45) と表される．ただし，ε_ijkは交代記号である．したがって，式(4.45)よりC₃₃は

C₃₃=−R

S (4.46)

と表され，C₃₃を除いて主不変量I，II，IIIを構成し直すと次式のようになる．

I⁰ =P−R

S (4.47)

II⁰ =Q−P R

S (4.48)

III⁰= 1 (4.49)

ただし

P =C_αα (4.50)

Q= 1

2(C_ααC_ββ−C_αβC_βα)−C_3αC_α3 (4.51)

R =ε_ijαC_i1C_j2C_α3−1 (4.52)

S =ε_αβ3C_α1C_β2 (4.53)

である．以上のように修正された主不変量I⁰，II⁰，III⁰を用いて，Ogden体の第2種

Piola-Kirchhoff応力テンソルを求めると次式のようになる．

S = 2 µ

a_I∂I⁰

∂C +a_II∂II⁰

∂C +a_III∂III⁰

∂C

¶

+ 2p∂III

∂C (4.54)

なお上式では，付録で示したI⁰，II⁰，III⁰のCによる微分の計算式を用いる．

ここで，C₃₃ はI⁰，II⁰，III⁰ に関して独立した変数であるため，(∂I⁰)/(∂C₃₃) = 0， (∂II⁰)/(∂C₃₃) = 0，(∂III⁰)/(∂C₃₃) = 0 となる．したがって，平面応力条件S³³ = 0 からp= 0が導かれる．以上より，式(4.54)は次式のように書き直される．

S= 2

· a_I

µ∂I

∂C −C⁻¹ S

¶ +a_II

µ∂II

∂C −PC⁻¹ S

¶¸

(4.55) また，構成則テンソルDは次式のようになる．

D = 4 µ

b_I ∂²I⁰

∂C∂C +b_II ∂²II⁰

∂C∂C +c_{I I}∂I⁰

∂C ⊗ ∂I⁰

∂C +c_{II II}∂II⁰

∂C ⊗∂II⁰

∂C +c_{I II}∂I⁰

∂C ⊗ ∂II⁰

∂C

¶

(4.56) 上式では，次式のI⁰，II⁰，III⁰のC による微分の計算式を用いる．

∂I⁰

∂C = ∂I

∂C −C⁻¹ S ,

∂II⁰

∂C = ∂II

∂C −PC⁻¹ S

∂III⁰

∂C = 0

∂²I⁰

∂C∂C = ∂²I

∂C∂C

−

ÃC⁻¹⊗C⁻¹+^∂C_∂C⁻¹

S −

∂C∂S ⊗C⁻¹+C⁻¹⊗_∂C^∂S S²

!

∂²II⁰

∂C∂C = ∂²II

∂C∂C

−P

ÃC⁻¹⊗C⁻¹+^∂C_∂C⁻¹

S −

∂C∂S ⊗C⁻¹+C⁻¹⊗ _∂C^∂S S²

!

− µ∂P

∂C ⊗C⁻¹

S +C⁻¹ S ⊗ ∂P

∂C

¶

∂²III⁰

∂C∂C = 0 (4.57)

ここで

∂C_ij⁻¹

∂C_kl = 1

2(C_ik⁻¹C_jl⁻¹+C_il⁻¹C_jk⁻¹)

∂S

∂C_kl =ε_3ikε_3jlC_ij

∂P

∂C_kl =δ_kl−δ_3kδ_3l (4.58)

である．ただし，δ_klはKroneckerのdeltaである．なお，a_K，b_K，c_{K L}は，4.2.1節で導 いた未知係数と同様のものを用いる．

4.2.5 未知係数aK，bK，cK Lの取り扱い

4.2.1節より，未知係数a_K，b_K，c_{K L}は，Λ_iの条件によって4通りの表現があり，非 圧縮条件式(4.6)まで考慮すると，Λ_iが2重根，3重根のときは一般にa_K6=b_Kであるこ とを導いた．一方，W を主不変量I，II，IIIの関数として，未知係数a_K，b_K，c_{K L}を 考慮した場合，次の関係式が恒等的に成り立つ．

a_I= ∂W

∂I , a_II = ∂W

∂II, a_III = ∂W

∂III, b_I = ∂W

∂I , b_II = ∂W

∂II, b_III = ∂W

∂III, c_{I I} = ∂²W

∂I∂I, c_{II II} = ∂²W

∂II∂II, c_{III III} = ∂²W

∂III∂III, c_{I II} = ∂²W

∂I∂II, c_{I III} = ∂²W

∂I∂III, c_{II III} = ∂²W

∂II∂III

(4.59)

式(4.59)より，明らかにa_K =b_Kが成り立つ．

一般に，Ogden体の構成則モデルは，非圧縮性まで考慮した場合，構成則の形は一意に

はならない．(例えば，式(4.6)より，ひずみエネルギー関数W の表現は一意にはならな い．) その結果，ひずみエネルギー関数W の微分後の形が異なる場合がある．文献^(4.12) においても，非圧縮超弾性体の構成則で用いる非圧縮の拘束関数の選択によって，応力テ ンソル・構成則テンソルの表現形式が変わり，また解の精度に差が現れることが報告され ている．応力テンソル・構成則テンソルを求める際に現れる未知係数a_K，b_Kも，恒等式 (4.59)によってa_K =b_Kが得られる一方，文献^(4.7)の方法に従って非圧縮条件下で求め た場合，a_K 6=b_Kの結果が得られるという，非圧縮性から生じる表現形式の多様性がみら れる．したがって，本研究において，未知係数a_K，b_K，c_{K L}の表現形式によって，解に どのような影響を及ぼすかを検証する．

ドキュメント内 目 次 (ページ 65-74)