Microsoft PowerPoint - 流体構造連成解析.ppt

流体–構造連成解析

東京大学大学院

新領域創成科学研究科

渡邉浩志

渡邉浩志 自己紹介

なぜ有限要素法を学ぶか

•

歴史の流れは、当初能力のある限られた一部の人のみが可能だったこ

とが、科学技術の発展により、一般人でも可能になる日が訪れる。という

方向にあると思われる。

•

身近になった有限要素法 ＣＡＤに組み込まれるものもある。

•

設計する ＝ 力学的に合理的である必要がある。

•

その意味で、設計と解析は不可分

•

設計担当者は外注のみで自分では解析せず、という状況でもないらしい

（簡単な解析は自分で行う）

•

汎用コード 解析データを与えさえすれば何らかの答えが出てくるが、こ

れが正しいかどうか判断するのはユーザである。

•

得てして重大な間違いに気がつかない場合が多い。

•

最近

FEMの品質保証

という言葉がよく聞かれるが、結局、やっていること

は、

設計・解析担当者の品質保証

にすぎないと思われる面がある。

有限要素法の正しい使い方

• 解析結果の間違いに 気がつかない ⇔ 気がつく

• 力学の経験に基づく皮膚感覚ともいえる、現象・モノに対す

る理解が必要

• 解析結果を正しく見る力があれば FEM は知らなくともよい

• ワープロが怪しい漢字変換をしても、日本語を知っているか

ら修正できる。

• これからエンジニアになる人にとっては、試作品をたくさん

作って壊す、という時代ではない。実験の機会も限られてい

る。（最後に一回だけなど。。。）

• FEMで経験を積んで、センスを磨くことには合理性がある。

• FEM の中身に対する理解があれば、より正しい使い方がわ

かる．また、より進んだ使い方がわかる．

• 現象・モノに対する謙虚な姿勢は必要（実験を伴わないシ

ミュレーションは無意味）

• 解析する能力とプログラミングの能力は原則として別物

• しかしながら、その成長には互いの協調が不可欠．

ヒトとサルの違い

• 道具を使う

• 道具とは？

• 代表が棒

中世と近代科学の差

（の一つ）

は材料力学

• 材料力学の創始者は

ガリレオ（二つの新科

学対話）

• ゴシック様式の建造は、

作っては壊れたの連

続だった

• 材料力学などの物理

学は数学者の興味を

引いた（Cauchy,

Euler, Lagrange,

Gauss・・・）

現在では？

• エンジニアなら知っていて当たり前？

• メシの種（by ７８７のチーフエンジニア）

• 計算機との関係は（有限要素法（冬学期））

• 材料力学の限界は、軽量化のための穴

• 連続体力学 → 美しい式（＝手計算では解けない）

• どっちが大切？

東京タワーの設計も、スカイツリーの

設計も設計も基本はトラスと梁

スカイツリーも最初は1本の梁で

解析した

• ただし断面は場所で変化する

• ちなみに、巨大タンカーも、飛行機も

思想は同じ

本日の講義内容

• 流体–構造連成解析の概要

• 連成手法の分類

– 連成手法の概要

– 一体型解法(強連成)

– 時差解法(弱連成)

– 分離型反復解法(漸近的強連成)

• 構造の有限要素法解析

– 境界値問題と変分原理

– 構成式

– 非線形方程式の解法

– 離散化

• ALE 有限要素法

– ALE 流体解析

– ALE 流体構造連成解析

• 解析例

流体–構造連成解析

• 流体–構造連成問題と

はたとえば紙が舞い落

ちる現象, 旗やこいの

ぼりが風にはためく現

象, あるいは心臓が拍

動し血液が流入出する

といった現象など,

流体

力が構造の変形をもた

らし, 同時に変位・変形

する構造が流れ場に影

響を及ぼす

ことを考慮

しなれけばならない力

学問題のことである..

流体–構造連成解析の概要

• 構造解析：Lagrange 表記、流体解析：Euler 表記、こ

れをどのように整合させるか

• 流体側で

境界面追跡型

の手法(interface-tracking

method) をとるか

境界面捕捉型

の手法(interface-capturing method) をとるか

• 流体と構造の方程式を解く段階で,

一体型解法

(monolithic method) と

分離型解法

(partitioned

method) のいずれをとるか

• 境界面での連成は

弱連成法

(weak coupling) と

強連

成法

_{(strong coupling)のいずれをとるか}

• 個々の問題に対して計算効率, 要求精度, 安定性の

観点から選択することになると考えられる.

境界面追跡型と境界面捕捉型

• ALE 法に代表される境界面追跡型の方法は流体

の移動境界を扱うため,

厳密に境界面が表現できる

反面, 構造が大きく変形した場合, 流体領域を離散

化した

_{メッシュが破綻し, 解析ができなくなる}

場合が

ある.

• これに対しImmersed Boundary Method,

Immersed Finite Element Method, Fictitious

Domain Methodなど境界面捕捉型の方法では

境

界面近傍の近似精度は劣り, 収束性や安定性が保

障されているわけではない

_{が, メッシュの破綻を考}

一体型解法

• 一体型解法とは流体と構造の連成系全体の方程式

を完全に連立させて強連成解法として解く手法のこ

とであり, 狭義の強連成法とも呼ばれる.

• 流体と構造の相互作用が強い場合, 他の方法に較

べて高い安定性・収束性を示す

• 解くべき方程式の元数が増え, 専用のプログラムを

開発する必要がある.

• また, 流体部は差分法, 構造部は有限要素法という

ように別々の離散化手法を採用することはできない.

分離型解法

• 分離型解法においては, 流体と構造の平衡方程式を個別に

立てて解く

• 何らかの反復計算によって両者の連成面での境界条件を整

合させる強連成法

• あるいは境界条件を厳密に満足させることなく時間進行させ

る弱連成法を採用するかの選択肢がある.

• 既存の流体解析コードと構造解析コードを活用することがで

き, 流体部に差分法を採用することも可能である. また並列

化も実装しやすい. しかしながら, 強い非線形性や連成効果

がある問題に対しては数値的に不安定になる場合もあるし,

原理的に適応できない現象もある.

連成問題の簡単な例題

• ２つの領域X, Y がそれぞれ１つの状態変数

を持ち, 以下のような時間発展型の

微分方程式で支配されるとする

強連成法と弱連成

• 抽象的な問題

• 強連成ではこれをそのまま解く

• 弱連成では本来未知数である

を

で置き換える

• 二つの式は独立に解くことができる

一体型解法

• [一体型解法]

を与えて, 次の連立方程式

を満たす

_を解く.

• 簡単な例題だと、下式をそのまま解くことに相当

• 連成効果の強い問題も適応可能。しかし大規模

な問題を解く必要があり、専用プログラムも必要

並列時差解法(弱連成)

• [並列時差解法]

を与えて, 次の２つの独立な方

程式を満たす

を解く

ここで,

はn ステップ以前のデータから陽的に

求められる

_{の予測値である.}

• 簡単な例題だと、たとえば

として

逐次時差解法

• [逐次時差解法]

1.

を与えて, 次の方程式を満たす

を解く.

ここで,

はn ステップ以前のデータから陽的に

求められる

の予測値

2.

を与えて, 次の方程式を満たす

を解く.

• 簡単な例題だと

分離型反復解法(漸近的強連成)

• 時差解法は弱連成法であるため連成効果を厳密に

近似できていないため,

強い非線形性や連成効果

を有する問題

においては

数値に不安定

になることが

知られている.

• 一体型解法も２つの問題を同様に有限要素法で離

散化するような場合はアルゴリズムの開発は比較

的容易であるが,問題の特性に応じて, たとえば有

限要素法と差分法などを組み合わせる必要が生じ

た場合には,

アルゴリズムの構築は困難

_となる.

• そこで, 近年では時差解法のような

分離型解法を基

本

_{とし, １つの計算ステップ内で}

反復計算

を行うこと

により, 強連成としての連成効果を漸近的満たす

分

離型反復解法

_{(iterative partitioned method) が注}

block Jacobi 法

• [block Jacobi 法] 以下をサブサイクルk に対して,

所定の条件を満たすまで反復する.

を与えて, 次の２つの独立な方程式

を満たす

_を解く.

• 簡単な例題だと

とし

て

構造解析

• 大学院講義資料から抜粋

• http://www.sml.k.u-tokyo.ac.jp/

members/nabe/lecture2007

あるいは

http://www.sml.k.u-tokyo.ac.jp/

members/nabe/lecture2005/

• 資料のほかにサンプルプログラムなども掲載

境界値問題

下図に示すような_, 物体 A についての境界値問題 _[B] を考える_. [B] 物体 A が占める領域を _{Ω, Ω} の境界を ∂Ω とし_{, ∂Ω} の部分集合 ∂Ω_D 上では変位境界条 件が与えられているものとする_. このような系に表面力 t, 体積力 ρg が作用するとき_, つり合い 条件を満たす変位 u ∈ V を求めよ_. ただし_{, ρ} は密度_{, g} は重力加速度_{, V} は変位の許容関数すな わち変位境界条件を満たす解の候補全体の集合とする_. 物体 A 表面力t 体積力 ρg ∂Ω_D ∂Ω 図 9: 境界値問題概念図 40

境界値問題の定式化

この問題は以下のように定式化できる

_.

[B]

与えられた

t, g

に対し

_,

以下を満たすような

u ∈ V

を求めよ

_.

[1 ]

平衡方程式

_(Cauchy

の第

₁

運動法則

₎

∇

x

· T + ρg = 0

(98)

∇

X

·



S · F

T



+ ρ

0

g = 0

(99)

[2 ]

境界条件式

u = u on ∂Ω

D

(100)

T

T

_{· n = t on ∂Ω − ∂Ω}

D

(101)



S · F

T



T

_{· N = ˜t}

(102)

[3 ]

変位・ひずみ関係式

[4 ]

応力・ひずみ関係式（構成式）

•

このうち

_{[1], [2]}

はいかなる問題でも共通（ただし，必要に応じて等価なもの

に書き換えることがある．

），

_[4]

は対象とする物質モデルで決まる．

_[3]

は

_[4]

に対応するものが選ばれる．

41

境界問題の弱形式定式化

1

• 平衡方程式と，力学的境界条件式を満足する一価連続な応力 T 全体の集合を T , 変位境界条 件を満足する一価連続な u 全体の集合をU とする_. • ˇT ∈ T , ˇu ∈ U とすると_, 下式が成立する_.  v (∇x · ˇT + ρg) · ˇu dv = 0 (103) • 上式に発散定理を用いると_, 以下のようになる_. ただし_{, s}t, su は t, u が指定されている s の 領域を示し_{, s = s}t + su である.  v ˇ T : (ˇu ⊗ ∇x) dv =  s_t t · ˇu ds +  s_u n · ˇT · u ds +  v ρg · ˇudv ₍₁₀₄₎ • ここで_, 一価連続で su 上で w = 0 となる w 全体の集合をW とする. • このとき u ∈ U, ∀w ∈ W_ˇ に対して u + w ∈ U_ˇ である_. • 平衡状態にある_Cauchy 応力を T として_, 下式が成立する_.  vT : (ˇu ⊗ ∇x) dv =  st t · ˇu ds +  su n · T · u ds +  v ρg · ˇudv (105)  vT : {(ˇu + w) ⊗ ∇x} dv =  st t · (ˇu + w) ds +  su n · T · u ds +  v ρg · (ˇu + w)dv (106) 42

境界問題の弱形式定式化

2

• 平衡状態にある_Cauchy 応力を T として_, 下式が成立する_.  vT : (ˇu ⊗ ∇x) dv =  st t · ˇu ds +  su n · T · u ds +  v ρg · ˇudv (107)  vT : {(ˇu + w) ⊗ ∇x} dv =  st t · (ˇu + w) ds +  su n · T · u ds +  v ρg · (ˇu + w)dv (108) • 両式を両辺それぞれ引くことにより下式が得られる_.  v T : (w ⊗ ∇x) dv =  s_t t · w ds +  v ρg · wdv ₍₁₀₉₎ • これは以下のように表すことができる_.  v T : δA(L) dv =  δv t · w ds +  v ρg · w dv ₍₁₁₀₎ ただし δA_(L) は_{, w ∈ W} から導かれる _Almange ひずみの線形成分である_. δA_{(L)ij =} 1 2  ∂wi ∂xj + ∂wj ∂xi (111) • また，この式の参照配置を変更することにより，下式が得られる．  V S : δE dv =  δV t · w dS +  V ρg · w dV ₍₁₁₂₎ 43

境界問題の弱形式定式化

3

 v T : δALdv  V S : δEdV (113)  V dV =  v 1 Jdv (114) T = 1 JF · S · F T _{→ S} ij = JF_im−1TmnF_jn−1 (115) δE = FTδA_LF (FT_δA LF )ij = FkiδALklFlj = ∂xk ∂X_i  1 2  ∂δu_k ∂x_l + ∂δu_l ∂x_k  ∂x_l ∂X_j = 1₂  ∂x_k ∂X_i ∂δu_k ∂x_l ∂x_l ∂X_j + ∂x_k ∂X_i ∂δu_l ∂x_k ∂x_l ∂X_j = 1₂  ∂x_k ∂X_i ∂δu_k ∂X_j + ∂x_l ∂X_i ∂δu_l ∂X_j = 1 2  δ_ki∂uk ∂X_i ∂δu_k ∂X_j + ∂δu_l ∂X_j  δ_li∂uk ∂X_i  = 1 2  ∂δu_i ∂X_j + ∂δu_j ∂X_i + ∂δu_k ∂X_i ∂u_k ∂X_j + ∂u_k ∂X_i ∂δu_k ∂X_j = δEij (116) 44

 V S : δEdV =  v (FT_δA LF ) : (JF−1 · T · F−T)_J1dv =  v J(F_kiδA_LklF_lj)(F_im−1T_mnF_jn−1)1 Jdv =  v δ_kmδ_lnA_LklT_mndv =  v T : δALdv (117) 45

total Lagrange

と

updated Lagrange 1

• 境界値問題は弱形式にすると  v T : δA(L) dv =  ∂v t · w ds +  v ρg · w dv (303)  V S : δE dV =  ∂V t · w dS +  V ρg· w dV (304) • 速度型にすると， _ v δA : ˙S_t(t) + 1 2 * δF_t(t)T · L + LT · δF_t(t) + : T dv = δ ˙R (305)  Ω S_ijδE_ijdV · =  Ω ˙S_ijδE_ij + S_ijδ ˙E_ijdΩ (306) • updated のところで使った構成則 ˙S_{t(t)ij = Cijkl}D_kl (307) • ˙St(t) は Truesdell 応力速度，相対 Kirchoff 応力のOldroyd 速度で，客観性がある．

• Total のところで使った構成式

S_ij = C_ijklE_kl (308) • Cijkl が一定とすると ( ˙Sij, ˙Ekl ともに，観測不変量)

˙S_{ij = Cijkl} ˙E_{kl (309)}

total Lagrange

と

updated Lagrange 2

• 異なる構成式 ˙ St(t)_ij = ¯CijklDkl (310) Sij = CijklEkl (311) ˙ Sij = CijklE˙kl (312) • 両者の関係は， ˙S0(t) = J0(t)F0(t)−1 ˙St(t)F0(t)−T (313) ˙ E0(t) = F0(t)TDF0(t) (314) を用いることにより， ¯ Cpqrs = 1 JFpiFqjFrkFslCijkl (315) と導かれる． • このような構成式の変換を行って初めて，２つの定式化は一致する 95

構成式

1

•

ここまでの理論は

_,

出発点の「連続体仮説」と力学の原理以外に特別な仮説を

導入することなく数学的に展開されてきた

_.

•

これに対し

_,

物質の力学的応答を記述する式

_,

すなわち構成式

_{(constitutive}

equa-tion)

についての理論は実験結果に基づく物理的考察から何らかの仮説を導入

せざるを得ない現象論的な側面を持つ

_.

•

しかしながら，実験式といっても数学的に正しくなくてなならない．

•

何らかの仮定をおいて組み立てた，数学的に正しいモデルが実験結果をよく

あわしているものが，一般に使われている構成式

34

構成式

2

•

変形した物体のある物質点における変形勾配が

F , Cauchy

応力が

T

であった

とする

_.

•

この物質点がさらに

Q

だけ剛体回転した場合

_,

テンソルの変換則により応力

は

Q · T · Q

T

と与えられる

_.

また

_,

このときの合計の変形勾配は

Q · F

である

_.

•

簡単のため

_{, T}

が

F

の関数として

T = f(F )

(78)

で与えられると仮定する

_.

•

ここで

f

はテンソルを変数とし

_,

関数値自身もテンソルとなるテンソル値テ

ンソル関数

_{(tensor–valued tensor function)}

である

_.

Q · T · Q

T

= f (Q · F )

(79)

あるいは

Q · f(F ) · Q

T

= f (Q · F )

(80)

が成立しなければならない

_.

このことから

f

の形は大きく限定される

_.

35

構成式

3

•

一般に

_,

構成式は次の

₃

つの原理に従わなければならない

_.

1.

応力決定の原理

_{(principle of determinism for the stress) :}

物体中の応力は物

体の運動履歴によって決定される

_.

2.

局所作用の原理

_{(principle of local action) :}

物質点

X

の応力の決定において

は

_{, X}

の近傍の運動だけが関与し

_,

その外の物質点の運動は無視できる

_.

3.

物質客観性の原理

_{(principle of material frame indiffernce or principle of}

mate-rial objectivity) :

構成式は基準枠

_{(referance frame)}

1

の変化の下で不変でなけ

ればならない

_.

すなわち

_{, 2}

つの基準枠

_O

∗

と

_O

が

x

∗

= c(t) + Q(t) · x

(81)

の関係にあるとき

T

∗

= Q(t) · T · Q(t)

T

(82)

の関係があり

_,

また構成式は両枠に対して同一の形をとらねばならない

_.

1_{事象(event) と呼ばれる空間内の位置 x と時間 t の組 {x, t} を観測する観測者 (observer) のこと. 基準枠が異なると事象は {x, t}, {x}∗_{, t}∗_{} とベクトル自体が異} なって観測される. 36

構成式

4

• 2つの基準枠 _O∗ と _O が x∗ = c(t) + Q(t) · x (83) • 2 点を結ぶベクトルは b∗ _{= x}∗₂ − x∗_{1, b = x2} − x₁ とすると， x∗₂ − x∗₁ = Q(t)(x2 − x1) (84) b∗ = Q(t) · b (85) このように両基準枠から観測されるベクトルが，限時刻での Q(t) のみを介した形で関係づ けられるとき，そのベクトルには客観性があるという． • テンソルはベクトルの線形変換なので，以下の関係を満たすとき客観性があるという K∗ = Q(t) · K · Q(t)T (86) • 客観性のあるテンソルの例 V , B, A, D, T (87) • 客観性のないテンソルの例 R, L, W (88) • 観測不変テンソルの例 U, C, E, S (89) 37

構成式

4

• 不合理な例：線形弾性体 εij(u) = 1 2  ∂ui ∂Xj + ∂uj ∂Xi (90) T_{ij = κ(div u)δij + 2Gε}D_{ij(u) (91)}

• x1 − x2 平面内で 90◦ 回転，応力は０のはず Fij = ⎡ ⎢ ⎣ 0 −1 0 1 0 0 0 0 1 ⎤ ⎥ ⎦ (92) Eij = 1 2 ⎛ ⎜ ⎝ ⎡ ⎢ ⎣ −1 −1 0 1 −1 0 0 0 0 ⎤ ⎥ ⎦ + ⎡ ⎢ ⎣ −1 1 0 −1 −1 0 0 0 0 ⎤ ⎥ ⎦ ⎞ ⎟ ⎠ = ⎡ ⎢ ⎣ −1 0 0 0 −1 0 0 0 0 ⎤ ⎥ ⎦ (93) • 微小変形の場合は，対角項は _{cos θ ≈ 1} 非対角項は_{sin θ ≈ θ} なので，問題ないが，有限変形 を考えると，不合理を生じる． 38

構成式

4

•

以下のような構成式を考えたとする．

T = aE

(94)

•

観測者が変わると

T

∗

_{= a}

∗

E

∗

となるが，

– T

∗

= QT Q

T

（客観性がある）

– a

∗

= a

（スカラー）

– E

∗

= E

（観測不変量）

•

従って

T

∗

= a

∗

E

∗

→ QT Q

T

= aE = T

(95)

となり，矛盾する．

•

たとえば，以下のものは矛盾しない．

T = aA

(96)

S = aE

(97)

39

長方形断面の梁要素

•

一般的な偏微分方程式を素直に離散化するなら、ソリッド要素が自然

•

しかし「構造解析」をする場合には、別な考え方がある

•

梁要素は代表的な構造要素

•

断面の寸法が長手方向の寸法より十分に小さいスレンダーな部材からなる構

造物の解析に使用される

_.

•

同構造をソリッド要素でモデル化するのに対して

_,

大幅に総自由度やモデル化

の手間が軽減される

_.

•

低次のソリッド要素で梁をモデル化すると，ロッキングする

•

時刻

₀

において，梁の中立軸に垂直な平面は変形の間も平面を保つが必ずし

も変形した中立軸に垂直である必要はない

•

梁の断面は変形しない．

4

座標及び変位の補間関数

1

• 時刻tで節点nの断面を張るベクトルを tV(n)₂ , tV(n)₃ , また時刻 t での節点 n の位置ベクトルを tx とする 全体座標系 局所座標系 a b e1 e2 e3 X₁ X₂ X₃ r1 r₂ r₃ t_V 2 t_V 3 1
1
2
2
3
3
4
4
図 1: 座標系の定義 • 時刻 t における要素内の任意の点における位置ベクトルは, t_x(r 1, r2, r3) = N(n)(r1)
t_xn₊ a 2r2tV (n) 2 + b 2r3tV (n) 3  (1) と, 表すことができる. • なお, 形状関数 N(n) には1次元要素のものを用いる. 5

座標及び変位の補間関数

2

• 局所座標ベクトル tV (n)_{2 ,} tV (n)₃ は梁の断面に固定されており_, 要素の変形とともに回転するこ とを考慮すれば_, 時刻 t から t_{(= t + ∆t)} までの間の変位増分ベクトル u は u = t
_{x −}t_x = N(n)  t
_x(n)₋t_x(n)₊ a 2r2  t
_V(n) 2 −tV(n)2  + b 2r3  t
_V(n) 3 −tV(n)3  となる_. • ここで∆t間の回転増分は微小であると仮定すれば_, t_V (n) i  t_V (n) i + θ(n) × t_V (n) i (2) と表すことができる_. • よって_, t
_V(n) i −tV (n) i = θ(n)×tV (n) i = ⎡ ⎢ ⎣ 0 −θ3(n) θ2(n) θ₃(n) 0 −θ(n)₁ −θ(n)₂ θ(n)₁ 0 ⎤ ⎥ ⎦ ⎧ ⎪ ⎨ ⎪ ⎩ t_V(n) i1 t_V(n) i2 t_V(n) i3 ⎫ ⎪ ⎬ ⎪ ⎭ = ⎧ ⎪ ⎨ ⎪ ⎩ θ₂(n)tV_i3(n) − θ₃(n)tV_i2(n) θ₃(n)tV_i1(n) − θ₁(n)tV_i3(n) θ₁(n)tV_i2(n) − θ₂(n)tV_i1(n) ⎫ ⎪ ⎬ ⎪ ⎭ = ⎡ ⎢ ⎣ 0 t_V(n) i3 −tV (n) i2 −t_V(n) i3 0 tV (n) i1 t_V(n) i2 −tV (n) i1 0 ⎤ ⎥ ⎦ ⎧ ⎪ ⎨ ⎪ ⎩ θ(n)₁ θ(n)₂ θ(n)₃ ⎫ ⎪ ⎬ ⎪ ⎭ (3) 6

座標及び変位の補間関数

3

• ここで tx, u の全体座標系での成分表示をしておく_. ⎧ ⎪ ⎨ ⎪ ⎩ t_x1 t_x2 t_x 3 ⎫ ⎪ ⎬ ⎪ ⎭ = N (n)_(r 1) ⎧ ⎪ ⎨ ⎪ ⎩ t_x(n) 1 t_x(n) 2 t_x(n) 3 ⎫ ⎪ ⎬ ⎪ ⎭+ a 2r2N(n)(r1) ⎧ ⎪ ⎨ ⎪ ⎩ t_V(n) 21 t_V(n) 22 t_V(n) 23 ⎫ ⎪ ⎬ ⎪ ⎭+ b 2r3N(n)(r1) ⎧ ⎪ ⎨ ⎪ ⎩ t_V(n) 31 t_V(n) 32 t_V(n) 33 ⎫ ⎪ ⎬ ⎪ ⎭ (4) ⎧ ⎪ ⎨ ⎪ ⎩ u₁ u₂ u₃ ⎫ ⎪ ⎬ ⎪ ⎭ = N (n)_(r 1) ⎧ ⎪ ⎨ ⎪ ⎩ u(n)₁ u(n)₂ u(n)₃ ⎫ ⎪ ⎬ ⎪ ⎭+ a 2r2N(n)(r1) ⎡ ⎢ ⎣ 0 t_V(n) 23 −tV22(n) −t_V(n) 23 0 tV21(n) t_V(n) 22 −tV21(n) 0 ⎤ ⎥ ⎦ ⎧ ⎪ ⎨ ⎪ ⎩ θ₁(n) θ₂(n) θ₃(n) ⎫ ⎪ ⎬ ⎪ ⎭ + b 2r3N(n)(r1) ⎡ ⎢ ⎣ 0 t_V(n) 33 −tV32(n) −t_V(n) 33 0 tV31(n) t_V(n) 32 −tV31(n) 0 ⎤ ⎥ ⎦ ⎧ ⎪ ⎨ ⎪ ⎩ θ₁(n) θ₂(n) θ₃(n) ⎫ ⎪ ⎬ ⎪ ⎭ (5) 7

シェル要素の定式化

• シェル要素の開発は過去多くの研究者により行われてきたが_, 現在も高精度_, 高信頼性かつ計 算効率のよい要素を目指して開発が続けられている_. 優れたシェル要素の条件は_, – 薄肉_, 厚肉双方のシェルに適用できる ₍薄肉シェル要素においてもロッキングが生じない_). – 任意の形状のシェルを取り扱える_. – 剛体モード以外の虚偽のゼロ固有値を持たない_. – 要素のゆがみに対して解の精度が損なわれない_. – 要素内の任意の点においてゼロを含む一定ひずみを表現できる_.

などが挙げられる_. これらを考慮すると_{, Bathe}ら_[?]が ₁₉₈₄ 年以降に発表した _{MITC (Mixed}

Interpolation of Tensorial Components) 要素が理論の明解さ_, 定式化の容易さの点で最も優れて

いる_.

幾何学的非線形シェル要素の基礎式

• Bathe らの _MITC 要素の概要を示す _{[?] [?].} この要素の特徴は_, 1. アイソパラメトリック _degenerated シェル要素である_. 2. 面外せん断ひずみに関しては_, あるサンプリング点の面外せん断ひずみの値から内挿する ように再定義する_. 3. ひずみ成分および応力成分には_, 自然座標系での共変成分_, 反変成分を使用する_.

などが挙げられる_. 上記の _{2, 3} が_{Mixed Interpolation of Tensorial Components} 要素と呼ばれる

理由であり_{, Assumed-Strain} 要素と参照されている場合もある_.

形状および変位補間式

• 図に時刻 ₀ におけるシェル要素の形状を示す_. 以下_, 特に記載のない場合には図中の記号で 説明をおこなう_. 図 ₅ より時刻 ₀ における要素内の任意の点の位置ベクトル X は形状関数 N (r1, r2₎ を用いて_, X = 4  n=1 Nnr1, r2Xn ₊ a 2r 34 n=1 Nnr1, r20V n_{3 (160)} で表される_. 式₍₁₆₀₎で 0V k₃ はシェルの厚み方向を示すベクトル₍以下_, ディレクター_. 左肩 の添字は時刻を表す₎であり_, 解析面に対する法線として一節点にただひとつ定義され_, その 節点を有する要素間で共有される_. ここでシェルの変形に関して次の ₂ つの仮定を設ける_. 1. 各節点におけるディレクター₍時刻 ₀ において直線を仮定した₎は_, 変形の間も直線を保つ が_, 必ずしも中立面に垂直である必要はない₍面外せん断変形を許す_). 2. シェルの肉厚 a は変形の間も変化しない₍大ひずみ領域は取り扱わない_). このとき時刻 t における要素内の任意の点 tx の位置ベクトルは_, t_{x =} 4  n=1 Nnr1, r2txn ₊ a 2r 34 n=1 Nnr1, r2tV n_{3 (161)} 52

となる_. したがって_, 時刻 t における変位 tu , 時刻 t から t_{(= t + ∆t)} までの変位増分 u は それぞれ tu = tx − X, u = tx −tx であるから_, t_{u =} 4  n=1 Nn_(r1, r2₎tUn ₊ a 2r 34 n=1 Nn_(r1, r2₎₍tV n₃ − 0V n_{3) (162)} u = 4  n=1 Nn_(r1, r2_)Un ₊ a 2r 34 n=1 Nn_(r1, r2₎₍tV n₃ − tV n_{3) (163)} となる_. ここで_, 各節点のシェルディレクターの時刻 t から t までの有限回転を表すテンソル を t_tR とすると_, tV k_{3 =} t_tR tV k₃ であるから式₍₁₆₃₎はつぎのように書き直される_. u = 4  n=1 Nn_(r1, r2_)Un ₊ a 2r 34 n=1 Nn_(r1, r2₎₍t_tR − I) tV n_{3 (164)} (I は単位テンソル₎ 53

a A B C D 1 2 3 4 e1 e2 e3 G1 G2 G3 x1 x2 x3 r1 r2 r3 α4 β4 u1 1 u1 2 u1 3 0_V4 1 0_V4 2 0_V4 3 a shell thickness xi Cartesian coordinates

ei orthonormal base vector in global coordinate system ri natural coordinates ( -1 ≤ ri ≤ 1)

Gi covariant base vector of natural coordinate ri

Vk

i orthonormal base vector in local coordinate system at node k

Vk

3 shell director vector

uk _{displacement vector at node k}

αk, βk rotations of the director vector about Vk₁ and Vk₂

図 5: 4 節点アイソパラメトリックシェル要素 54

面外せん断ひずみ成分補間式

• 通常のアイソパラメトリック要素では_, 要素内任意点のひずみ_, ひずみ増分は式_{(162), (163)}で 表せる変位_, 変位増分の空間微分をとることにより直接求められるが_{, MITC} 要素では_, 面外 せん断ひずみの成分についてのみ_, あるサンプリング点における面外せん断ひずみ ₍通常の方 法から求める₎ から内挿関数を再定義して求める_. 図 ₆ にサンプリング点位置 ₍点_A∼_D) を_, 式(165)に面外せん断ひずみ成分の補間式を示す. E_{13 =} 1 2(1 + r 2_)EA 13 + 1₂(1 − r2)E13C (165) E_{23 =} 1 2(1 + r 1_)ED 23 + 1 2(1 − r 1_)EB 23 (166) ここで_, E₁₃, E_{23 :} 要素内の任意点のせん断ひずみ EA∼D _: サンプリング点におけるせん断ひずみ である. 55

1 2 3 4 A A B B C C D D r1 r1 r1 r1 r1 r2 r2 r2 r 2 r2 EA 13 EC 13 EB 23 ED 23 図 6: 面外せん断ひずみのサンプリング点 56

ALE 流体解析

• Lagrange 表示で表される座標系を物質座標系

• Euler 表示で表される座標系を空間座標系

• これらの座標系とは独立した任意の座標系，すなわ

ち参照座標系を設定し，その座標系上で物体の運

動を記述することを参照表示(ALE 表示)と定義する．

• 今，Lagrange，Euler，ALE の各表示により解析領

域および領域内のある点を次のように表す

時間導関数の関係

• 実質時間導関数と物質時間導関数の関係式

• 実質時間導関数(物質時間導関数) と空間時間導関数の

関係式

• 実質時間導関数物(物質時間導関数) と参照時間導関数

の関係式

各種の速度の物理的解釈

• Euler 座標系に対する物質点の速度v．

• 参照座標系に対する物質点の速度w．

Euler 領域での時間導関数の式

• 位置ベクトルをx とすれば，速度関係式は

• 参照座標系に対する物質点の相対速度

を導入

すると

• これより

参考までに

Euler だと

ALE 流体構造連成解析（１）

• Navier-Stokes 方程式と連続の式のALE 表示

• 非圧縮性Newton 流体（構成則）

ALE 流体構造連成解析（２）

• 連成面での力学的平衡条件式と幾何学的連

続条件式

有限要素離散化したマトリックス方程式

• 流体

連成系のマトリックス方程式

• ただし肩符号c (common)は連成面上の自由

度で、肩符号i (independent)はそれ以外の

自由度とする

補足：大まかなプ

ログラムの流れ

• メッシュの制御

• 時間積分

補足：メッシュ制御

• 流体部分を構造が取り囲んでいるような系、total submerged system の場合、連成系のマトリックスを解く と構造の変位が得られる。 • 流体領域を弾性体に見立てて、流体と固体の境界上の 節点は、固体の解析で得られた変位を、変位境界条件と して与えて、弾性体の剛性方程式を解き、内部の節点の 移動量を計算する。 • 流体の節点のメッシュの速度は、単純に変位をΔt で割 るだけでもよい。（経験的にはそれほど差は無い） • なお、弾性体でなく、Laplace 方程式でも良いし、線形に 配分しても良いし、もっと複雑な構成式（超弾性体）などを 使っても良い。 • ただし、逆に言うと、 ALEはそれで対応できる範囲に限ら れると考えたほうがよい。構造の変形にともなって。流体 メッシュが大きくゆがんだり、極めて間隔が狭くなるような 場合は対応できない。メッシュを切りなおして、再補間す るなどの処理が必要。

？

補足：時間積分は？

• 固体を含まない流体の有限要素法では、速

度と加速度の項だけ。

• 固体があると、速度、加速度に加えて、変位

が未知数になるが、どうやって積分するのか。

• 一般的には陰解法で Newmark-β法をつか

う。

動的解析手法

• 静的な仮想仕事の式 _ V S_ijδE_ijdV =  V g_iρδu_idV +  S t_iδu_idS (1) 加速度の項がある場合 _ V ρa_iδu_idV +  V S_ijδE_ijdV =  V g_iρδu_idV +  S t_iδu_idS (2) • 時刻t+ Δtでの運動方程式 M ·t+Δt_{¨u + C ·}t+Δt_{˙u +}t+Δt_{Q =} t+Δt_{F (3)} • M は質量マトリックス，C は減衰マトリックス. 実はC は仮想仕事の式からは導かれない。現実問題では、固 定部の摩擦などによって減衰が生じる。また材料も粘性を持つのでこれからも減衰を生じる。これをあわせて C としているので実際には条件を変えて実験を行って C を妥当な値に調整する、ということが行われている。 • ここでは，この方程式をそのままの形で時間方向に積分する直接時間積分法についてのべる． • 線形問題では，固有モードの重ね合わせで解析する場合もあり，モード重ね合わせ法という． • 時刻 t での解が得られているとき，この解を基にして t+ Δt の解を厳密に求める方法を陰解法（線形加速度法， Newmark-β 法など），t の解を基に予測する方法を陽解法（中央差分）という • 通常は M, C は時刻によらず一定．しかし Q は変化する • 内力ベクトルについて線形化し，Newton-Raphson 法を用いる．

M ·t+Δt_¨u(k) _{+ C ·}t+Δt_˙u(k)₊t+Δt_K(k−1)_Δu(k) ₌ t+Δt_{F −}t+Δt_Q(k−1) ₍₄₎

t+Δt_u(k) ₌ t_u(k−1) _{+ Δu}(k) ₍₅₎

t+Δt_˙u(k) ₌ t _˙u(k−1) _{+ Δ ˙u}(k) ₍₆₎

t+Δt_¨u(k) ₌ t_¨u(k−1) _{+ Δ¨u}(k) ₍₇₎

• このままでは，未知量として Δu(k), Δ ˙u(k), Δ¨u(k) が含まれるため，通常は，この三者を何らかの仮定の下に関 係づけ，どれか１つの変数にまとめて解析する．

線形加速度法

• 線形加速度法は実際の解析では用いられることは少ない • 文字通り時刻 t から _{t + Δt} までの加速度が線形に変化すると仮定する． t+τ ¨ u = t ¨ u + τ Δt( t+τ ¨ u − t ¨ u) (8) • この式を τ を変数として積分し _{τ = Δt} を代入すると t+Δt ˙ u = t ˙ u + Δt 2 _t+Δt ¨ u +t ¨ u (9) t+Δt_{u =} t_{u + Δt} t ˙ u + Δt2 3 t ¨ u + Δt2 6 t+Δt ¨ u (10) • このとき， Δ ˙u(k) = Δt 2 Δ ¨u (k) ₍₁₁₎ Δu(k) = Δt 2 6 Δ ¨u (k) ₍₁₂₎ 6

Newmark-β

法

離散化された動的仮想仕事式から M · t ¨ u + C · t ˙ u + t_{Q =} t_F (13) 内力ベクトルtQを線形化し，_{Newton-Raphson}法を用いる． M · t ¨ u(k) + C ·tu˙ (k) +tK(k−1)Δu(k) = tF − tQ(k−1) (14) t_u(k) = tu(k−1) + Δu(k) (15) t ˙ u(k) = tu˙ (k−1) + Δ ˙u(k) (16) t ¨ u(k) = tu¨(k−1) + Δ ¨u(k) (17) k = 0では，_{t = 0}における収束値． t_K(0) = tK (18) t_Q(0) = tQ (19) t_U(0) = tU (20) t_U˙ (0) = tU˙ (21) t_U¨ (0) = tU¨ (22) (14)に₍₁₅₎∼₍₁₇₎を代入する． M · (t ¨ u(k−1) + Δ ¨u(k)) + C · (tu˙ (k−1) + Δ ˙u(k)) + tK(k−1)Δu(k) = tF − tQ(k−1) (23) 7

未知項を左辺に，既知項を右辺に移項

M · Δ¨u(k) + C · Δ ˙u(k) +tK(k−1)Δu(k) = tF − tQ(k−1) − M ·tu¨(k−1) − C ·tu˙ (k−1) (24) (24)で_{k = 1}のときは

M · Δ¨u(1) + C · Δ ˙u(1) +tKΔu(1) = tF − tQ − M · tu − C ·¨ tu˙ (25) Newmark-β 法の関係式 t ˙ u = t ˙ u + Δtγ tu + (1 − γ)_¨ tu_¨  (26) t_{u =} t_{u + Δt} t ˙ u + Δt2  β tu +_¨  1 2 − β t ¨ u (27) プライマリー変数を加速度にした場合． 1回目のイタレーションでは_(26),(27)を直接変形して，_{Δ ˙u}(1)_{, Δu}(1)を_{Δ ¨}u(1)であらわす Δ ˙u(1) = tu −˙ tu = Δt˙  γ tu + (1 − γ)_¨ tu_¨  (28) = Δt  γ(tu −¨ tu) +¨ tu¨  (29) = γΔtΔ ¨u(1) + Δt tu¨ (30) Δu(1) = tu − tu = Δt tu + Δt˙ 2  β tu +_¨  1 2 − β t ¨ u (31) = Δtu + Δt˙ 2  β  t ¨ u − t ¨ u + 1 2 t ¨ u (32) = βΔt2Δ ¨u(1) + Δt tu +˙ 1 2Δt 2 t_{u¨ (33)} 8

2回目以降のイタレーションではtu_˙ (k)とtu_˙ (k−1)をそれぞれ₍₂₆₎を用いてあらわし，差を取って Δ ˙u(k)を_{Δ ¨}u(k)であらわす．_Δu(k)についても同様． t ˙ u(k) = tu + Δt˙  γ tu_¨(k) _{+ (1 − γ)}tu_¨  (34) t ˙ u(k−1) = tu + Δt˙  γ tu_¨(k−1) _{+ (1 − γ)}tu_¨  (35) Δ ˙u(k) = tu˙ (k) − tu˙ (k−1) = γΔt  t ¨ u(k) − t ¨ u(k−1) (36) = γΔtΔ ¨u(k) (37) t_u(k) = tu + Δttu + Δt˙ 2  β tu_¨(k) ₊  1 2 − β t ¨ u (38) t_u(k−1) = tu + Δttu + Δt˙ 2  β tu_¨(k−1) ₊  1 2 − β t ¨ u (39) Δu(k) = tu(k) − tu(k−1) = βΔt2  t ¨ u(k) −t ¨ u(k−1) (40) = βΔt2Δ ¨u(k) (41) 以上により求められた_{Δ ˙u}(k)_{, Δu}(k)を離散化された増分形の運動方程式に代入する_. 1回目のイタレーション (25)に_(30),(33)を代入する

M · Δ¨u(1) + C · Δ ˙u(1) +tKΔu(1) = tF −tQ − M · tu − C ·¨ tu˙ (13)

M·Δ¨u(1)+C·(γΔtΔ ¨u(1)+Δttu)+¨ tK·(βΔt2Δ ¨u(1)+Δttu+˙ 1 2Δt 2t_{u) =¨} t_{F −}t_Q−M·t ¨ u−C·t ˙ u (42) 9

M · Δ¨u(1) + γΔtC · Δ ¨u(1) + βΔt2 tK · Δ¨u(1) = tF − tQ − M · tu − C ·¨ tu˙ − ΔtC ·t ¨ u − Δt t_{K ·} t ˙ u − 1 2Δt 2 t_{K ·}t ¨ u (43)  M + γΔtC + βΔt2 tKΔ ¨u(1) = tF −tQ − M · tu − C ·¨ tu˙ − ΔtC · t ¨ u − Δt t_{K ·} t ˙ u − 1 2Δt 2 t_{K ·} t ¨ u (44) 2回目以降のイタレーション (24)に_(37),(41)を代入する

M · Δ¨u(k) + C · Δ ˙u(k) +tK(k−1)Δu(k) = tF − tQ(k−1) − M ·tu¨(k−1) − C ·tu˙ (k−1) (12)

M ·Δ¨u(k)+C · (γΔtΔ ¨u(k)) +tK(k−1)·(βΔt2Δ ¨u(k)) = tF −tQ(k−1)−M ·tu¨(k−1)−C ·tu˙ (k−1) (45) (M + γΔtC + βΔt2 tK(k−1))Δ ¨u(k) = tF −tQ(k−1) − M · tu¨(k−1) − C · tu˙ (k−1) (46)

プライマリー変数を変位にする．

1回目のイタレーション

前述のように_Newmark-β法の関係式から，_{Δ ˙u}(1)_{, Δu}(1)を_{Δ ¨}u(1)であらわすと式_(30),(33)のように

なる Δ ˙u(1) = γΔtΔ ¨u(1) + Δt tu¨ (18) Δu(1) = βΔt2Δ ¨u(1) + Δttu +˙ 1 2Δt 2 t_{u¨ (21)} (33)から_{Δ ¨}u(1)を_Δu(1)で表す， βΔt2Δ ¨u(1) = Δu(1) − Δt tu −˙ 1 2Δt 2 t_{u¨ (47)} 10

Δ ¨u(1) = 1 βΔt2Δu (1) ₋ Δt βΔt2 t ˙ u − Δt2 2βΔt2 t ¨ u (48) = 1 βΔt2Δu (1) ₋ 1 βΔt t ˙ u − 1 2β t ¨ u (49) (49)を₍₃₀₎に代入して，_{Δ ˙u}(1)を_Δu(1)で表す， Δ ˙u(1) = γΔt  1 βΔt2Δu (1) ₋ 1 βΔt t ˙ u − 1 2β t ¨ u + Δt tu¨ (50) = γΔt βΔt2Δu (1) ₋ γΔt βΔt t ˙ u − γΔt 2β t ¨ u + Δt t ¨ u (51) = γ βΔtΔu (1) ₋ γ β t ˙ u − γΔt 2β t ¨ u + Δtt ¨ u (52) 当然ながら，_(26),(27)から直接計算することもできる． (27)から_{Δ ¨}u(1)を_Δu(1)で表す． t_{u =} t_{u + Δt} t ˙ u + Δt2  β tu +_¨  1 2 − β t ¨ u (15) t_{u −} t_{u = Δt}t ˙ u + Δt2  β  t ¨ u − t ¨ u + 1 2 t ¨ u (53) Δu(1) = Δt tu + Δt˙ 2βΔ ¨u(1) + Δt 2 2 t ¨ u (54) Δ ¨u(1) = 1 βΔt2 Δu (1) ₋ 1 βΔt t ˙ u − 1 2β t ¨ u (55) (26)から_{Δ ˙u}(1)を_{Δ ¨}u(1)で表したものが₍₃₀₎である． Δ ˙u(1) = γΔtΔ ¨u(1) + Δt tu¨ (18) 11

(30)に₍₅₅₎を代入する． Δ ˙u(1) = γΔt  1 βΔt2 Δu (1) ₋ 1 βΔt t ˙ u − 1 2β t ¨ u + Δt tu¨ (56) = γ βΔt Δu (1) ₋ γ β t ˙ u − γΔt 2β t ¨ u + Δtt ¨ u (57) (49)と₍₅₅₎，₍₅₂₎と₍₅₇₎は等しいことが確認できる． 2回目以降のイタレーション． Δ ˙u(k) = γΔtΔ ¨u(k) (25) Δu(k) = βΔt2Δ ¨u(k) (29) これより Δ ¨u(k) = 1 βΔt2 Δu (k) ₍₅₈₎ Δ ˙u(k) = γΔt βΔt2 Δu (k) ₌ γ βΔt Δu (k) ₍₅₉₎ 以上により求められた_{Δ ¨}u(1)_{, Δ ˙u}(1)を離散化された増分形の運動方程式に代入する． 1回目のイタレーション (25)に_(49),(52)を代入する

M · Δ¨u(1) + C · Δ ˙u(1) +tKΔu(1) = tF − tQ − M · tu − C ·¨ tu˙ (13)

M ·  1 βΔt2Δu (1) ₋ 1 βΔt t ˙ u − 1 2β t ¨ u + C ·  γ βΔtΔu (1) ₋ γ β t ˙ u − γΔt 2β t ¨ u + Δt t ¨ u +tKΔu(1) = tF − tQ − M · tu − C ·¨ tu˙ (60) 12

1 βΔt2M · Δu (1) ₊ γ βΔtC · Δu (1) ₊t_KΔu(1) = tF −tQ − M · tu − C ·¨ tu˙ + 1 βΔtM · t ˙ u + 1 2βM · t ¨ u + γ βC · t ˙ u + γΔt 2β C · t ¨ u − ΔtC · t ¨ u (61)  1 βΔt2M + γ βΔtC + t_K Δu(1) (62) = tF − tQ − M · tu − C ·¨ tu +˙ 1 βΔtM · t ˙ u + 1 2βM · t ¨ u + γ βC · t ˙ u + γΔt 2β C · t ¨ u − ΔtC · t ¨ u (63) = tF − tQ − M · tu − C ·¨ tu +˙  1 2βM +  γΔt 2β − Δt C t ¨ u +  1 βΔtM + γ βC t ˙ u (64) 確認のため加速度をプライマリー変数とする式₍₄₄₎に，_{Δ ¨}u(1)を_Δu(1)で表した式₍₄₉₎を代入して プライマリー変数を変位にする．  M + γΔtC + βΔt2 tKΔ ¨u(1) = tF −tQ − M · tu − C ·¨ tu − ΔtC ·˙ tu − Δt¨ tK · tu −˙ 1 2Δt 2 t_{K ·} t ¨ u (32)  M + γΔtC + βΔt2 tK  1 βΔt2Δu (1) ₋ 1 βΔt t ˙ u − 1 2β t ¨ u = tF −tQ − M · tu − C ·¨ tu − ΔtC ·˙ tu − Δt¨ tK · tu −˙ 1 2Δt 2 t_{K ·} t ¨ u (65) 13

 1 βΔt2M + γ βΔtC + t_K Δu(1) (66) = tF − tQ − M · tu − C ·¨ tu − ΔtC ·˙ tu − Δt¨ tK · tu −˙ 1 2Δt 2 t_{K ·} t ¨ u + 1 βΔtM · t ˙ u + γ βC t ˙ u + Δt tK · tu +˙ 1 2βM · t ¨ u + γΔt 2β C · t ¨ u + Δt2 2 t_{K ·} t ¨ u (67) = tF − tQ − M · tu − C ·¨ tu − ΔtC ·˙ tu +¨ 1 βΔtM · t ˙ u + γ βC t ˙ u + 1 2βM · t ¨ u + γΔt 2β C · t ¨ u (68) = tF − tQ − M · tu − C ·¨ tu +˙  1 2βM +  γΔt 2β − Δt C t ¨ u +  1 βΔtM + γ βC t ˙ u (69) (69)と₍₆₄₎は一致することが確認できる． 2回目以降のイタレーション． (24)に(58),(59)を代入する

M · Δ¨u(k) + C · Δ ˙u(k) +tK(k−1)Δu(k) = tF − tQ(k−1) − M ·tu¨(k−1) − C ·tu˙ (k−1) (12)  1 βΔt2M + γ βΔtC + t_K(k−1) Δu(k) = tF − tQ(k−1) − M · tu¨(k−1) − C ·tu˙ (k−1) (70) 確認のため₍₄₆₎に₍₅₈₎を代入する (M + γΔtC + βΔt2 tK(k−1))Δ ¨u(k) = tF −tQ(k−1) − M · tu¨(k−1) − C ·tu˙ (k−1) (34)  1 βΔt2M + γ βΔtC + t_K(k−1) Δu(k) = tF − tQ(k−1) − M · tu¨(k−1) − C ·tu˙ (k−1) (71) 14

更なる補足：非圧縮性拘束条件付の

超弾性体の場合

⎥

⎦

⎤

⎢

⎣

⎡

=

0

0

0

M

M

_⎥

⎦

⎤

⎢

⎣

⎡

=

0

T

B

B

K

K

•

質量マトリックスは、変位の項に値があるだけで、Lagrange 未定乗数のところ

は０

•

解けるのか？

•

Cの形にもよるが、よく用いられる C＝αM＋βKの形だと、係数が掛かるだ

けで、静的解析の場合と同じ式を解くことになる。

解析例ー心臓シミュレータ

•

心臓はモデル化に際してマルチフィジックスの現象を取り扱う必要がある挑戦的

なテーマであり、これまで様々な観点から多くの研究がなされてきた

•

分子・細胞科学が明らかにするミクロレベルの知見から、一人の人間という個体レベル

の間の相互作用を予測することは不可能

•

個々の実験結果を統合し一つのシステムとして計算機内に再現する

•

著者らは既報において、心筋の興奮伝播、興奮収縮連関を基礎とした左心室の流

体構造連成解析を行った （Watanabe, H. et. al, Biophys. J., vol. 87 . (2004)）

マルチスケール

マ

ル

チ

フ

ィ

ジ

ッ

ク

ス

分子構造・機能

細胞

臓器

生体

力学的現象

電気的現象

生化学反応

イオンチャンネル 活動電位 心電図 収縮蛋白 心筋張力 圧ー容積関係 K+ Na+ _Ca++ ブドウ糖 解糖 ATP 乳酸