流 体 構 造 連 成 解 析 東 京 大 学 大 学 院 新 領 域 創 成 科 学 研 究 科 渡 邉 浩 志
渡 邉 浩 志 自 己 紹 介 固 体 構 造 解 析 から 流 体 構 造 連 成 解 析 へ
なぜ 有 限 要 素 法 を 学 ぶか 歴 史 の 流 れは 当 初 能 力 のある 限 られた 一 部 の 人 のみが 可 能 だったこ とが 科 学 技 術 の 発 展 により 一 般 人 でも 可 能 になる 日 が 訪 れる という 方 向 にあると 思 われる 身 近 になった 有 限 要 素 法 CADに 組 み 込 まれるものもある 設 計 する = 力 学 的 に 合 理 的 である 必 要 がある その 意 味 で 設 計 と 解 析 は 不 可 分 設 計 担 当 者 は 外 注 のみで 自 分 では 解 析 せず という 状 況 でもないらしい ( 簡 単 な 解 析 は 自 分 で 行 う) 汎 用 コード 解 析 データを 与 えさえすれば 何 らかの 答 えが 出 てくるが こ れが 正 しいかどうか 判 断 するのはユーザである 得 てして 重 大 な 間 違 いに 気 がつかない 場 合 が 多 い 最 近 FEMの 品 質 保 証 という 言 葉 がよく 聞 かれるが 結 局 やっていること は 設 計 解 析 担 当 者 の 品 質 保 証 にすぎないと 思 われる 面 がある
有 限 要 素 法 の 正 しい 使 い 方 解 析 結 果 の 間 違 いに 気 がつかない 気 がつく 力 学 の 経 験 に 基 づく 皮 膚 感 覚 ともいえる 現 象 モノに 対 す る 理 解 が 必 要 解 析 結 果 を 正 しく 見 る 力 があれば FEM は 知 らなくともよい ワープロが 怪 しい 漢 字 変 換 をしても 日 本 語 を 知 っているか ら 修 正 できる これからエンジニアになる 人 にとっては 試 作 品 をたくさん 作 って 壊 す という 時 代 ではない 実 験 の 機 会 も 限 られてい る ( 最 後 に 一 回 だけなど ) FEMで 経 験 を 積 んで センスを 磨 くことには 合 理 性 がある FEM の 中 身 に 対 する 理 解 があれば より 正 しい 使 い 方 がわ かる.また より 進 んだ 使 い 方 がわかる. 現 象 モノに 対 する 謙 虚 な 姿 勢 は 必 要 ( 実 験 を 伴 わないシ ミュレーションは 無 意 味 ) 解 析 する 能 力 とプログラミングの 能 力 は 原 則 として 別 物 しかしながら その 成 長 には 互 いの 協 調 が 不 可 欠.
ヒトとサルの 違 い 道 具 を 使 う 道 具 とは? 代 表 が 棒 棒 は 引 張 圧 縮 曲 げどれが 強 いか
中 世 と 近 代 科 学 の 差 (の 一 つ)は 材 料 力 学 材 料 力 学 の 創 始 者 は ガリレオ( 二 つの 新 科 学 対 話 ) ゴシック 様 式 の 建 造 は 作 っては 壊 れたの 連 続 だった 材 料 力 学 などの 物 理 学 は 数 学 者 の 興 味 を 引 いた(Cauchy, Euler, Lagrange, Gauss )
現 在 では? エンジニアなら 知 っていて 当 たり 前? メシの 種 (by 787のチーフエンジニア) 計 算 機 との 関 係 は( 有 限 要 素 法 ( 冬 学 期 )) 材 料 力 学 の 限 界 は 軽 量 化 のための 穴 連 続 体 力 学 美 しい 式 (= 手 計 算 では 解 けない) どっちが 大 切?
東 京 タワーの 設 計 も スカイツリーの 設 計 も 設 計 も 基 本 はトラスと 梁
スカイツリーの 仕 口 の 解 析
スカイツリーも 最 初 は1 本 の 梁 で 解 析 した ただし 断 面 は 場 所 で 変 化 する ちなみに 巨 大 タンカーも 飛 行 機 も 思 想 は 同 じ
本 日 の 講 義 内 容 流 体 構 造 連 成 解 析 の 概 要 連 成 手 法 の 分 類 連 成 手 法 の 概 要 一 体 型 解 法 ( 強 連 成 ) 時 差 解 法 ( 弱 連 成 ) 分 離 型 反 復 解 法 ( 漸 近 的 強 連 成 ) 構 造 の 有 限 要 素 法 解 析 境 界 値 問 題 と 変 分 原 理 構 成 式 非 線 形 方 程 式 の 解 法 離 散 化 ALE 有 限 要 素 法 ALE 流 体 解 析 ALE 流 体 構 造 連 成 解 析 解 析 例
流 体 構 造 連 成 解 析 流 体 構 造 連 成 問 題 と はたとえば 紙 が 舞 い 落 ちる 現 象, 旗 やこいの ぼりが 風 にはためく 現 象, あるいは 心 臓 が 拍 動 し 血 液 が 流 入 出 する といった 現 象 など, 流 体 力 が 構 造 の 変 形 をもた らし, 同 時 に 変 位 変 形 する 構 造 が 流 れ 場 に 影 響 を 及 ぼすことを 考 慮 しなれけばならない 力 学 問 題 のことである..
流 体 構 造 連 成 解 析 の 概 要 構 造 解 析 :Lagrange 表 記 流 体 解 析 :Euler 表 記 こ れをどのように 整 合 させるか 流 体 側 で 境 界 面 追 跡 型 の 手 法 (interface-tracking method) をとるか 境 界 面 捕 捉 型 の 手 法 (interfacecapturing method) をとるか 流 体 と 構 造 の 方 程 式 を 解 く 段 階 で, 一 体 型 解 法 (monolithic method) と 分 離 型 解 法 (partitioned method) のいずれをとるか 境 界 面 での 連 成 は 弱 連 成 法 (weak coupling) と 強 連 成 法 (strong coupling)のいずれをとるか 個 々の 問 題 に 対 して 計 算 効 率, 要 求 精 度, 安 定 性 の 観 点 から 選 択 することになると 考 えられる.
境 界 面 追 跡 型 と 境 界 面 捕 捉 型 ALE 法 に 代 表 される 境 界 面 追 跡 型 の 方 法 は 流 体 の 移 動 境 界 を 扱 うため, 厳 密 に 境 界 面 が 表 現 できる 反 面, 構 造 が 大 きく 変 形 した 場 合, 流 体 領 域 を 離 散 化 したメッシュが 破 綻 し, 解 析 ができなくなる 場 合 が ある. これに 対 しImmersed Boundary Method, Immersed Finite Element Method, Fictitious Domain Methodなど 境 界 面 捕 捉 型 の 方 法 では 境 界 面 近 傍 の 近 似 精 度 は 劣 り, 収 束 性 や 安 定 性 が 保 障 されているわけではないが, メッシュの 破 綻 を 考 慮 する 必 要 のない 方 法 を 構 成 できる.
一 体 型 解 法 一 体 型 解 法 とは 流 体 と 構 造 の 連 成 系 全 体 の 方 程 式 を 完 全 に 連 立 させて 強 連 成 解 法 として 解 く 手 法 のこ とであり, 狭 義 の 強 連 成 法 とも 呼 ばれる. 流 体 と 構 造 の 相 互 作 用 が 強 い 場 合, 他 の 方 法 に 較 べて 高 い 安 定 性 収 束 性 を 示 す 解 くべき 方 程 式 の 元 数 が 増 え, 専 用 のプログラムを 開 発 する 必 要 がある. また, 流 体 部 は 差 分 法, 構 造 部 は 有 限 要 素 法 という ように 別 々の 離 散 化 手 法 を 採 用 することはできない.
分 離 型 解 法 分 離 型 解 法 においては, 流 体 と 構 造 の 平 衡 方 程 式 を 個 別 に 立 てて 解 く 何 らかの 反 復 計 算 によって 両 者 の 連 成 面 での 境 界 条 件 を 整 合 させる 強 連 成 法 あるいは 境 界 条 件 を 厳 密 に 満 足 させることなく 時 間 進 行 させ る 弱 連 成 法 を 採 用 するかの 選 択 肢 がある. 既 存 の 流 体 解 析 コードと 構 造 解 析 コードを 活 用 することがで き, 流 体 部 に 差 分 法 を 採 用 することも 可 能 である. また 並 列 化 も 実 装 しやすい. しかしながら, 強 い 非 線 形 性 や 連 成 効 果 がある 問 題 に 対 しては 数 値 的 に 不 安 定 になる 場 合 もあるし, 原 理 的 に 適 応 できない 現 象 もある.
連 成 問 題 の 簡 単 な 例 題 2つの 領 域 X, Y がそれぞれ1つの 状 態 変 数 を 持 ち, 以 下 のような 時 間 発 展 型 の 微 分 方 程 式 で 支 配 されるとする 後 退 Eular 型 時 間 積 分 を 適 用 すると
強 連 成 法 と 弱 連 成 抽 象 的 な 問 題 強 連 成 ではこれをそのまま 解 く 弱 連 成 では 本 来 未 知 数 である で 置 き 換 える を 二 つの 式 は 独 立 に 解 くことができる
一 体 型 解 法 [ 一 体 型 解 法 ] を 与 えて, 次 の 連 立 方 程 式 を 満 たす を 解 く. 簡 単 な 例 題 だと 下 式 をそのまま 解 くことに 相 当 連 成 効 果 の 強 い 問 題 も 適 応 可 能 しかし 大 規 模 な 問 題 を 解 く 必 要 があり 専 用 プログラムも 必 要
並 列 時 差 解 法 ( 弱 連 成 ) [ 並 列 時 差 解 法 ] 程 式 を 満 たす を 与 えて, 次 の2つの 独 立 な 方 を 解 く ここで, はn ステップ 以 前 のデータから 陽 的 に 求 められる の 予 測 値 である. 簡 単 な 例 題 だと たとえば として
逐 次 時 差 解 法 [ 逐 次 時 差 解 法 ] 1. を 与 えて, 次 の 方 程 式 を 満 たす を 解 く. ここで, はn ステップ 以 前 のデータから 陽 的 に 求 められる の 予 測 値 2. を 与 えて, 次 の 方 程 式 を 満 たす を 解 く. 簡 単 な 例 題 だと
分 離 型 反 復 解 法 ( 漸 近 的 強 連 成 ) 時 差 解 法 は 弱 連 成 法 であるため 連 成 効 果 を 厳 密 に 近 似 できていないため, 強 い 非 線 形 性 や 連 成 効 果 を 有 する 問 題 においては 数 値 に 不 安 定 になることが 知 られている. 一 体 型 解 法 も2つの 問 題 を 同 様 に 有 限 要 素 法 で 離 散 化 するような 場 合 はアルゴリズムの 開 発 は 比 較 的 容 易 であるが, 問 題 の 特 性 に 応 じて, たとえば 有 限 要 素 法 と 差 分 法 などを 組 み 合 わせる 必 要 が 生 じ た 場 合 には, アルゴリズムの 構 築 は 困 難 となる. そこで, 近 年 では 時 差 解 法 のような 分 離 型 解 法 を 基 本 とし, 1つの 計 算 ステップ 内 で 反 復 計 算 を 行 うこと により, 強 連 成 としての 連 成 効 果 を 漸 近 的 満 たす 分 離 型 反 復 解 法 (iterative partitioned method) が 注 目 されている.
block Jacobi 法 [block Jacobi 法 ] 以 下 をサブサイクルk に 対 して, 所 定 の 条 件 を 満 たすまで 反 復 する. を 与 えて, 次 の2つの 独 立 な 方 程 式 を 満 たす を 解 く. 簡 単 な 例 題 だと て とし
構 造 解 析 大 学 院 講 義 資 料 から 抜 粋 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 Ω D ρg Ω 9: 40
. [B] t, g, u V. [1 ] (Cauchy 1 ) 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, ǔ U,. ( x Ť + ρg) ǔ dv = 0 (103) v,., s t,s u t, u s, s = s t + s u. Ť :(ǔ x )dv = t ǔ ds + n Ť u ds + v s t s u, s u w =0 w W. ǔ U, w W ǔ + w U. Cauchy T,. v v ρg ǔdv (104) T :(ǔ x )dv = t ǔ ds + n T u ds + ρg ǔdv (105) v s t s u v T : {(ǔ + w) x } dv = t (ǔ + w)ds + n T u ds + ρg (ǔ + w)dv (106) s t s u v 42
2 Cauchy T,. v T :(ǔ x )dv = t ǔ ds + n T u ds + ρg ǔdv (107) v s t s u v T : {(ǔ + w) x } dv = t (ǔ + w)ds + n T u ds + ρg (ǔ + w)dv (108) s t s u v. T :(w x )dv = t w ds + v s t. T : δa (L) dv = v δv t w ds + δa (L), w W Almange. δa (L)ij = 1 ( wi + w ) j 2 x j x i v v ρg wdv (109) ρg w dv (110) (111) S : δe dv = t w ds + ρg w dv (112) V δv V 43
3 v T : δa L dv V dv = V v S : δedv (113) 1 dv (114) J T = 1 J F S F T S ij = JF 1 im T mnf 1 jn (115) δe = F T δa L F (F T δa L F ) ij = F ki δa Lkl F lj = x { ( k 1 δuk + δu )} l xl X i 2 x l x k X j = 1 ( xk δu k x l + x ) k δu l x l 2 X i x l X j X i x k X j = 1 ( xk δu k + x ) l δu l 2 X i X j X i X j = 1 {( ) u k δuk δ ki + δu ( )} l u k δ li 2 X i X j X j X i = 1 ( δui + δu j + δu k u k + u ) k δu k 2 X j X i X i X j X i X j = δe ij (116) 44
V S : δedv = = = = v v v v (F T δa L F ):(JF 1 T F T ) 1 J dv J(F ki δa Lkl F lj )(F 1 im T mnf 1 jn ) 1 J dv δ km δ ln A Lkl T mn dv T : δa L dv (117) 45
total Lagrange updated Lagrange 1 updated v v V T : δa (L) dv = S : δe dv = δa : Ṡt(t)+ 1 2 ( Ω S ij δe ij dv V v t w ds + t w ds + V v ρg w dv (303) ρg w dv (304) ( ) δf t (t) T L + L T δf t (t) : T dv = δṙ (305) ) = Ω Ṡ ij δe ij + S ij δėijdω (306) Ṡ t (t) ij = C ijkl D kl (307) Ṡt(t) Truesdell Kirchoff Oldroyd Total S ij = C ijkl E kl (308) C ijkl (Ṡij, Ėkl ) Ṡ ij = C ijkl Ė kl (309) 94
total Lagrange updated Lagrange 2 Ṡ t (t) ij = C ijkl D kl (310) S ij = C ijkl E kl (311) Ṡ ij = C ijkl Ė kl (312) Ṡ 0 (t) =J 0 (t)f 0 (t) 1 Ṡ t (t)f 0 (t) T (313) Ė 0 (t) =F 0 (t) T DF 0 (t) (314) C pqrs = 1 J F pif qj F rk F sl C ijkl (315) 95
1,.,, (constitutive equation). 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, 3. 1. (principle of determinism for the stress) :. 2. (principle of local action) : X, X,. 3. (principle of material frame indiffernce or principle of material 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 2 x 1, b = x 2 x 1 x 2 x 1 = Q(t)(x 2 x 1 ) (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 ( ui + u ) j 2 X j X i (90) T ij = κ(div u)δ ij +2Gε D ij(u) (91) x 1 x 2 90 0 1 0 F ij = 1 0 0 (92) 0 0 1 E ij = 1 1 1 0 1 1 0 1 1 0 + 1 1 0 2 0 0 0 0 0 0 1 0 0 = 0 1 0 (93) 0 0 0 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
.,. 0 4
1 t n t V (n) 2, t V (n) 3, t n t x t V 3 b 1 3 4 t V 2 2 r 3 r 2 a r X 1 3 1 3 4 2 X 2 e 2 e 3 e 1 X 1 1: t, ( t x(r 1,r 2,r 3 )=N (n) (r 1 ) t x n + a 2 r 2 t V (n) 2 + b ) 2 r 3 t V (n) 3 (1),., N (n) 1. 5
2 t V (n) 2, t V (n) 3,, t t (= t + t) u. u = t x t x = N (n) {t x (n) t x (n) + a 2 r 2 t, ( ) t V (n) 2 t V (n) 2 + b ( ) } 2 r t 3 V (n) 3 t V (n) 3., t V (n) i t V (n) i t V (n) i t V (n) i = θ (n) t V (n) i 0 θ (n) 3 θ (n) 2 t V (n) = θ (n) 3 0 θ (n) i1 1 θ (n) 2 θ (n) t V (n) i2 1 0 t V (n) = i3 0 t V (n) i3 t V (n) i2 θ (n) = t V (n) i3 0 t V (n) 1 i1 θ (n) t V (n) i2 t V (n) 2 i1 0 θ (n) 3 + θ (n) t V (n) i (2) θ (n) 2 t V (n) i3 θ (n) 3 t V (n) i1 θ (n) 1 t V (n) i2 θ (n) θ (n) 3 t V (n) i2 1 t V (n) i3 θ (n) 2 t V (n) i1 (3) 6
3 t x, u. t x 1 t x 2 t x 3 u 1 u 2 u 3 = N (n) (r 1 ) = N (n) (r 1 ) t x (n) 1 t x (n) 2 t x (n) 3 u (n) 1 u (n) 2 u (n) 3 + a t V (n) 21 2 r 2N (n) (r 1 ) t V (n) 22 t V (n) + b 2 r 3N (n) (r 1 ) 23 + a 0 t V (n) 23 t V (n) 22 2 r 2N (n) (r 1 ) t V (n) 23 0 t V (n) 21 t V (n) 22 t V (n) 21 0 + b 0 t V (n) 33 t V (n) 32 2 r 3N (n) (r 1 ) t V (n) 33 0 t V (n) 31 t V (n) 32 t V (n) 31 0 t V (n) 31 t V (n) 32 t V (n) 33 θ (n) 1 θ (n) 2 θ (n) 3 θ (n) 1 θ (n) 2 θ (n) 3 (4) (5) 7
,,.,, ( )......, Bathe [?] 1984 MITC (Mixed Interpolation of Tensorial Components),. 50
Bathe MITC [?] [?]., 1. degenerated. 2.,. 3.,,.. 2, 3 Mixed Interpolation of Tensorial Components, Assumed-Strain. 51
0.,. 5 0 X N(r 1,r 2 ), X = 4 N ( n r 1,r 2) X n + a 2 r3 n=1 4 N ( n r 1,r 2) 0 V n 3 (160) n=1. (160) 0 V k 3 (,. ),,. 2. 1. ( 0 ),, ( ). 2. a ( ). t t x, t x = 4 N ( n r 1,r 2) t x n + a 2 r3 n=1 4 N ( n r 1,r 2) t V n 3 (161) n=1 52
., t t u, t t (= t + t) u t u = t x X, u = t x t x, t u = u = 4 N n (r 1,r 2 ) t U n + a 2 r3 n=1 4 N n (r 1,r 2 )U n + a 2 r3 n=1 4 N n (r 1,r 2 )( t V n 3 0 V n 3) (162) n=1 4 N n (r 1,r 2 )( t V n 3 t V n 3 ) (163)., t t t t R, t V k 3 = t t R t V k 3 (163). 4 u = N n (r 1,r 2 )U n + a 4 2 r3 N n (r 1,r 2 )( t t R I) t V n 3 (164) n=1 n=1 n=1 (I ) 53
r 3 r 2 B 2 A u 1 3 3 G 3 G 2 G 1 u 1 1 1 u 1 2 a 0 V 4 3 D r 1 x 3 C 0 V 4 1 4 β 4 0 V 4 2 e 3 α 4 e 2 e 1 x 2 x 1 a shell thickness x i Cartesian coordinates e i orthonormal base vector in global coordinate system r i natural coordinates ( -1 r i 1) covariant base vector of natural coordinate r i G i V k i V k 3 u k orthonormal base vector in local coordinate system at node k shell director vector displacement vector at node k α k, β k rotations of the director vector about V k 1 and V k 2 5: 4 54
,, (162), (163),, MITC,, ( ). 6 ( A D), (165). E 13 = 1 2 (1 + r2 )E13 A + 1 2 (1 r2 )E13 C (165) E 23 = 1 2 (1 + r1 )E23 D + 1 2 (1 r1 )E23 B (166), E 13,E 23 : E A D :. 55
6: r 2 2 A 1 B D r 1 3 C 4 E A 13 r 2 r 2 A r 1 r 1 E C 13 C r 2 r 2 D E D 23 E B 23 r 1 r 1 B 56
ALE 流 体 解 析 Lagrange 表 示 で 表 される 座 標 系 を 物 質 座 標 系 Euler 表 示 で 表 される 座 標 系 を 空 間 座 標 系 これらの 座 標 系 とは 独 立 した 任 意 の 座 標 系,すなわ ち 参 照 座 標 系 を 設 定 し,その 座 標 系 上 で 物 体 の 運 動 を 記 述 することを 参 照 表 示 (ALE 表 示 )と 定 義 する. 今,Lagrange,Euler,ALE の 各 表 示 により 解 析 領 域 および 領 域 内 のある 点 を 次 のように 表 す
時 間 導 関 数 の 関 係 実 質 時 間 導 関 数 と 物 質 時 間 導 関 数 の 関 係 式 実 質 時 間 導 関 数 ( 物 質 時 間 導 関 数 ) と 空 間 時 間 導 関 数 の 関 係 式 実 質 時 間 導 関 数 物 ( 物 質 時 間 導 関 数 ) と 参 照 時 間 導 関 数 の 関 係 式
各 種 の 速 度 の 物 理 的 解 釈 Euler 座 標 系 に 対 する 物 質 点 の 速 度 v. 参 照 座 標 系 に 対 する 物 質 点 の 速 度 w. Euler 座 標 系 に 対 する 参 照 座 標 系 の 速 度
Euler 領 域 での 時 間 導 関 数 の 式 位 置 ベクトルをx とすれば, 速 度 関 係 式 は 参 照 座 標 系 に 対 する 物 質 点 の 相 対 速 度 すると これより を 導 入 参 考 までに Euler だと
ALE 流 体 構 造 連 成 解 析 (1) Navier-Stokes 方 程 式 と 連 続 の 式 のALE 表 示 非 圧 縮 性 Newton 流 体 ( 構 成 則 ) 弾 性 体 の 運 動 方 程 式
ALE 流 体 構 造 連 成 解 析 (2) 連 成 面 での 力 学 的 平 衡 条 件 式 と 幾 何 学 的 連 続 条 件 式 連 成 系 の 弱 形 式
有 限 要 素 離 散 化 したマトリックス 方 程 式 流 体 構 造
連 成 系 のマトリックス 方 程 式 ただし 肩 符 号 c (common)は 連 成 面 上 の 自 由 度 で 肩 符 号 i (independent)はそれ 以 外 の 自 由 度 とする
マトリックスの 成 分
補 足 : 大 まかなプ ログラムの 流 れ メッシュの 制 御 時 間 積 分
補 足 :メッシュ 制 御 流 体 部 分 を 構 造 が 取 り 囲 んでいるような 系 total submerged system の 場 合 連 成 系 のマトリックスを 解 く と 構 造 の 変 位 が 得 られる? 流 体 領 域 を 弾 性 体 に 見 立 てて 流 体 と 固 体 の 境 界 上 の 節 点 は 固 体 の 解 析 で 得 られた 変 位 を 変 位 境 界 条 件 と して 与 えて 弾 性 体 の 剛 性 方 程 式 を 解 き 内 部 の 節 点 の 移 動 量 を 計 算 する 流 体 の 節 点 のメッシュの 速 度 は 単 純 に 変 位 をΔt で 割 るだけでもよい ( 経 験 的 にはそれほど 差 は 無 い) なお 弾 性 体 でなく Laplace 方 程 式 でも 良 いし 線 形 に 配 分 しても 良 いし もっと 複 雑 な 構 成 式 ( 超 弾 性 体 )などを 使 っても 良 い ただし 逆 に 言 うと ALEはそれで 対 応 できる 範 囲 に 限 ら れると 考 えたほうがよい 構 造 の 変 形 にともなって 流 体 メッシュが 大 きくゆがんだり 極 めて 間 隔 が 狭 くなるような 場 合 は 対 応 できない メッシュを 切 りなおして 再 補 間 す るなどの 処 理 が 必 要
補 足 : 時 間 積 分 は? 固 体 を 含 まない 流 体 の 有 限 要 素 法 では 速 度 と 加 速 度 の 項 だけ 固 体 があると 速 度 加 速 度 に 加 えて 変 位 が 未 知 数 になるが どうやって 積 分 するのか 一 般 的 には 陰 解 法 で Newmark-β 法 をつか う
t +Δt V S ij δe ij dv = g i ρδu i dv + t i δu i ds (1) V V S ρa i δu i dv + S ij δe ij dv = g i ρδu i dv + t i δu i ds (2) V V M t+δt ü + C t+δt u + t+δt Q = t+δt F (3) S M C. C C C t t +Δt Newmark-β t M, C Q Newton-Raphson M t+δt ü (k) + C t+δt u (k) + t+δt K (k 1) Δu (k) = t+δt F t+δt Q (k 1) (4) t+δt u (k) = t u (k 1) +Δu (k) (5) t+δt u (k) = t u (k 1) +Δ u (k) (6) t+δtü (k) = t ü (k 1) +Δü (k) (7) 4
Δu (k),δ u (k),δü (k) 5
t t +Δt t+τü = t ü + τ Δt (t+τ ü t ü) (8) τ τ =Δt t+δt u = t u + Δt 2 t+δt u = t u +Δt t u + Δt2 3 ( t+δtü + t ü ) (9) tü + Δt2 t+δtü (10) 6 Δ u (k) = Δt 2 Δü(k) (11) Δu (k) = Δt2 6 Δü(k) (12) 6
Newmark-β t Q Newton-Raphson M t ü + C t u + t Q = t F (13) M t ü (k) + C t u (k) + t K (k 1) Δu (k) = t F t Q (k 1) (14) t u (k) = t u (k 1) +Δu (k) (15) t u (k) = t u (k 1) +Δ u (k) (16) t ü (k) = t ü (k 1) +Δü (k) (17) k =0 t =0 t K (0) = t K (18) t Q (0) = t Q (19) t U (0) = t U (20) t U (0) = t U (21) t Ü (0) = t Ü (22) (14) (15) (17) M ( t ü (k 1) +Δü (k) )+C ( t u (k 1) +Δ u (k) )+ t K (k 1) Δu (k) = t F t Q (k 1) (23) 7
M Δü (k) + C Δ u (k) + t K (k 1) Δu (k) = t F t Q (k 1) M t ü (k 1) C t u (k 1) (24) (24) k =1 M Δü (1) + C Δ u (1) + t KΔu (1) = t F t Q M tü C t u (25) Newmark-β { } t u = t u +Δt γ t ü +(1 γ) t ü { ( ) t 1 u = t u +Δt t u +Δt 2 β t ü + 2 β 1 (26),(27) Δ u (1), Δu (1) Δü (1) { } Δ u (1) = t u t u =Δt γ t ü +(1 γ) t ü { } =Δt γ( t ü t ü)+ t ü tü } (26) (27) (28) (29) = γδtδü (1) +Δt t ü (30) ( ) 1 Δu (1) = t u t u =Δt t u +Δt {β 2 t ü + 2 β tü } (31) ( ) =Δ t u +Δt {β 2 t ü t ü + 1 } (32) 2tü = βδt 2 Δü (1) +Δt t u + 1 2 Δt2 t ü (33) 8
2 t u (k) t u (k 1) (26) Δ u (k) Δü (k) Δu (k) { } t u (k) = t u +Δt γ t ü (k) +(1 γ) t ü { } t u (k 1) = t u +Δt γ t ü (k 1) +(1 γ) t ü ) (t Δ u (k) = t u (k) t u (k 1) = γδt ü (k) t ü (k 1) (34) (35) (36) = γδtδü (k) (37) { ( ) t 1 u (k) = t u +Δt t u +Δt 2 β t ü (k) + 2 β tü } (38) { ( ) t 1 u (k 1) = t u +Δt t u +Δt 2 β t ü (k 1) + 2 β tü } (39) ( ) Δu (k) = t u (k) t u (k 1) = βδt 2 t ü (k) t ü (k 1) (40) = βδt 2 Δü (k) (41) Δ u (k), Δu (k). 1 (25) (30),(33) M Δü (1) + C Δ u (1) + t KΔu (1) = t F t Q M tü C t u (13) M Δü (1) +C (γδtδü (1) +Δt t ü)+ t K (βδt 2 Δü (1) +Δt t u+ 1 2 Δt2t ü)= t F t Q M tü C t u (42) 9
M Δü (1) + γδtc Δü (1) + βδt 2 t K Δü (1) = t F t Q M tü C t u ΔtC tü Δt t K t u 1 2 Δt2 t K tü (43) ( M + γδtc + βδt 2 t K ) Δü (1) = t F t Q M tü C t u ΔtC tü Δt t K t u 1 2 Δt2 t K tü (44) 2 (24) (37),(41) M Δü (k) + C Δ u (k) + t K (k 1) Δu (k) = t F t Q (k 1) M t ü (k 1) C t u (k 1) (12) M Δü (k) +C (γδtδü (k) )+ t K (k 1) (βδt 2 Δü (k) )= t F t Q (k 1) M tü (k 1) C t u (k 1) (45) (M + γδtc + βδt 2 t K (k 1) )Δü (k) = t F t Q (k 1) M tü (k 1) C t u (k 1) (46) 1 Newmark-β Δ u (1), Δu (1) Δü (1) (30),(33) Δ u (1) = γδtδü (1) +Δt t ü (18) Δu (1) = βδt 2 Δü (1) +Δt t u + 1 2 Δt2 t ü (21) (33) Δü (1) Δu (1) βδt 2 Δü (1) =Δu (1) Δt t u 1 2 Δt2 t ü (47) 10
Δü (1) = 1 βδt 2Δu(1) = 1 βδt 2Δu(1) 1 βδt Δt βδt 2 t u t u 1 2β Δt2 2βΔt 2 (48) tü (49) (49) (30) Δ u (1) Δu (1) ( 1 Δ u (1) = γδt 1 u 1 ) +Δt t ü (50) βδt 2Δu(1) βδt 2β = γδt γδt u γδt +Δt t ü βδt 2Δu(1) βδt 2β (51) = γ βδt Δu(1) γ u γδt +Δt t ü β 2β (52) (26),(27) (27) Δü (1) Δu (1) { ( ) t 1 u = t u +Δt t u +Δt 2 β t ü + 2 β { ( ) t u t u =Δt t u +Δt 2 t β ü t ü + 1 tü } 2 Δu (1) =Δt t u +Δt 2 βδü (1) + Δt2 Δü (1) = 1 βδt 2 Δu(1) 1 βδt (26) Δ u (1) Δü (1) (30) 2 t u 1 2β } tü (15) (53) tü (54) tü (55) Δ u (1) = γδtδü (1) +Δt t ü (18) 11
(30) (55) Δ u (1) = γδt ( 1 βδt 2 Δu(1) 1 βδt = γ βδt Δu(1) γ β t u γδt 2β (49) (55) (52) (57) 2 t u 1 2β tü ) +Δt t ü (56) tü +Δt t ü (57) Δ u (k) = γδtδü (k) (25) Δu (k) = βδt 2 Δü (k) (29) Δü (k) = 1 βδt 2 Δu(k) (58) Δ u (k) = γδt βδt 2 Δu(k) = γ βδt Δu(k) (59) Δü (1), Δ u (1) 1 (25) (49),(52) M M Δü (1) + C Δ u (1) + t KΔu (1) = t F t Q M tü C t u (13) ( 1 1 t u 1 tü ) ( γ + C βδt 2Δu(1) βδt 2β βδt Δu(1) γ t u γδt tü ) +Δt t ü + t KΔu (1) β 2β = t F t Q M tü C t u 12 (60)
( 1 βδt 2M + 1 βδt 2M Δu(1) + γ βδt C Δu(1) + t KΔu (1) = t F t Q M tü C t u + 1 βδt M t u + 1 2β M tü + γ β C t u + γδt 2β C tü ΔtC tü (61) γ βδt C + t K ) Δu (1) = t F t Q M tü C t u + 1 βδt M t u + 1 2β M tü + γ β C t u + γδt 2β C tü ΔtC tü = t F t Q M tü C t u + { ( ) } 1 γδt 2β M + 2β Δt C tü ( 1 + βδt M + γ ) β C (62) (63) t u (64) (44) Δü (1) Δu (1) (49) ( M + γδtc + βδt 2 t K ) Δü (1) = t F t Q M tü C t u ΔtC tü Δt t K t u 1 2 Δt2 t K tü (32) ( M + γδtc + βδt 2 t K ) ( 1 βδt 2Δu(1) 1 βδt t u 1 2β tü ) = t F t Q M tü C t u ΔtC tü Δt t K t u 1 2 Δt2 t K tü (65) 13
( 1 βδt 2M + γ ) βδt C + t K Δu (1) = t F t Q M tü C t u ΔtC tü Δt t K t u 1 2 Δt2 t K tü + 1 βδt M t u + γ β C t u +Δt t K t u + 1 2β M tü + γδt 2β C tü + Δt2 t K tü (67) 2 = t F t Q M tü C t u ΔtC tü + 1 βδt M t u + γ β C t u + 1 2β M tü + γδt 2β C tü (68) = t F t Q M tü { ( ) } 1 γδt C t u + 2β M + 2β Δt C tü ( 1 + βδt M + γ ) β C t u (69) (69) (64) 2 (24) (58),(59) M Δü (k) + C Δ u (k) + t K (k 1) Δu (k) = t F t Q (k 1) M t ü (k 1) C t u (k 1) (12) ( 1 βδt 2M + γ ) βδt C + t K (k 1) Δu (k) = t F t Q (k 1) M t ü (k 1) C t u (k 1) (70) (46) (58) (M + γδtc + βδt 2 t K (k 1) )Δü (k) = t F t Q (k 1) M tü (k 1) C t u (k 1) (34) ( 1 βδt 2M + γ ) βδt C + t K (k 1) Δu (k) = t F t Q (k 1) M t ü (k 1) C t u (k 1) (71) (66) 14
更 なる 補 足 : 非 圧 縮 性 拘 束 条 件 付 の 超 弾 性 体 の 場 合 M = M 0 0 0 K = K B T B 0 質 量 マトリックスは 変 位 の 項 に 値 があるだけで Lagrange 未 定 乗 数 のところ は0 解 けるのか? Cの 形 にもよるが よく 用 いられる C=αM+βKの 形 だと 係 数 が 掛 かるだ けで 静 的 解 析 の 場 合 と 同 じ 式 を 解 くことになる
マ ル チ フ ィ ジ ッ ク ス 解 析 例 ー 心 臓 シミュレータ 力 学 的 現 象 生 化 学 反 応 電 気 的 現 象 ブドウ 糖 解 糖 乳 酸 ATP 収 縮 蛋 白 イオンチャンネル Na + Ca ++ 心 筋 張 力 活 動 電 位 分 子 構 造 機 能 細 胞 臓 器 生 体 マルチスケール 圧 ー 容 積 関 係 心 電 図 心 臓 はモデル 化 に 際 してマルチフィジックスの 現 象 を 取 り 扱 う 必 要 がある 挑 戦 的 なテーマであり これまで 様 々な 観 点 から 多 くの 研 究 がなされてきた 分 子 細 胞 科 学 が 明 らかにするミクロレベルの 知 見 から 一 人 の 人 間 という 個 体 レベル の 間 の 相 互 作 用 を 予 測 することは 不 可 能 個 々の 実 験 結 果 を 統 合 し 一 つのシステムとして 計 算 機 内 に 再 現 する 著 者 らは 既 報 において 心 筋 の 興 奮 伝 播 興 奮 収 縮 連 関 を 基 礎 とした 左 心 室 の 流 体 構 造 連 成 解 析 を 行 った (Watanabe, H. et. al, Biophys. J., vol. 87. (2004)) K +
aortic valve aorta left atrium 心 周 期 Mitral valve left ventricle pressure 収 縮 ( 拍 出 ) 末 期 Phase2 phase3 phase1 phase4 Volume 左 心 室 には 大 動 脈 弁 僧 房 弁 があり 自 己 収 縮 弛 緩 す る Phase 1 等 容 収 縮, (0.1 sec ) Phase 2 拍 出, (0.2 sec) ( 約 55 ~ 60% が 拍 出 される) Phase 3 等 容 弛 緩 ( 0.1 sec ) Phase 4 流 入 (0.6 sec) 拡 張 ( 流 入 ) 末 期
Peterson four state model シミュレータの 概 要 Caイオン 濃 度 と 収 縮 力 の 関 係 RV RP LP RAV LAV モノドメインモデル 興 奮 伝 播 心 臓 電 気 現 象 FHNモデル 細 胞 興 奮 に 関 する 最 も 簡 略 なモデル 興 奮 - 収 縮 関 連 心 臓 機 械 現 象 LinYinのモデル 流 体 構 造 連 成 心 筋 の 運 動 血 液 の 影 響 Modified Alexander model Windkessel Model PV Structure CP FRL P2 C R2 1 ELA P1 FAo R1 P 40 Aortic Valve Fluid FMi 10 Mitral Valve 要 素 数 : 心 筋 9792, 左 心 室 内 部 流 体 18976 ともにQ1P0 総 自 由 度 : 約 120000
FSI analysis Fluid:ALE form of N-S equation Automatic mesh update (solving elastic equation) Standard SUPG stabilization Structure:total Lagrange formulation Finite element formulation Fluid Structure : strong coupling method Fluid mesh : weak coupling method Time integration : Predictor multi corrector algorithm Solver : GMRES with ILU precondition
心 室 内 血 流 の 可 視 化 心 室 内 部 の 流 体 構 造 を 観 察 す るために 導 入 したパーティクル トレースでは 流 入 時 に 僧 帽 弁 に 相 当 する 面 の 流 量 に 比 例 す るようにパーティクルを 生 成 し ている 流 入 した 血 液 は 直 後 の 拍 出 期 には 拍 出 されず 代 わりに 古 い 血 液 が 拍 出 される 古 い 血 液 が 心 室 内 に 停 留 する ことなく 更 新 されていくことは, 血 栓 生 成 の 防 止 に 繋 がると 考 えられ 合 理 的 であるといえる.
渦 の 役 割 (1)
渦 の 役 割 (2)