FEBio：生体組織や生体材料の解析に特化した非線形有限要素解析ソフトウェア

By M.Sato

概要

流体と構造の連成解析を行う場合、一つのプログラムでこれらの解析を実行できれば非常に便利である。 このような機能を有するプログラムにFEBioというソフトウエアがある。

FEBioは、生体組織や生体材料の解析に特化した非線形有限要素解析ソフトウェアであり、 オープンソースとしてソースコードが公開されているFEBio。

今回はFEBioについて調査を行いテスト計算を行ってみたのでその結果を記述する。

FEBioの機能

FEBioは、生体組織や生体材料の挙動を高精度にシミュレーションするために設計された有限要素解析ソフトウェアで、 多様な解析モジュールと物理現象に対応している。

構造解析機能

  • 非線形構造解析
    大変形・大ひずみ・非線形材料特性を持つ構造の挙動を解析可能。生体組織のような高非線形挙動にも対応。

  • 弾性・超弾性材料モデル
    Neo-Hookean、Mooney-Rivlin、Ogden、Holzapfel-Gasser-Ogdenなど、生体材料に適した様々な弾性・超弾性材料モデルを搭載。

  • 塑性・損傷モデル
    弾塑性、粘塑性、損傷進展を含む材料挙動をシミュレーション可能。

  • 粘弾性材料特性
    時間依存的な材料挙動を再現。Prony系列などによるモデル化が可能で、軟部組織などの挙動に対応。

  • 異方性材料
    線維構造を持つ材料の方向依存性（例：筋肉や腱）を再現する異方性モデルを提供。

  • 接触解析
    摩擦あり/なしの接触定義、スライディング、分離、嵌合など、複雑な接触状況を精密に解析。

マルチフィジックス解析機能

  • 熱伝導解析
    生体内の温度分布や熱伝導挙動の解析が可能。

  • 拡散解析（化学種輸送）
    濃度勾配に基づく物質拡散（例：酸素や栄養分）をシミュレート。

  • 多相材料モデル（poroelasticity）
    固体と流体の両方の挙動を扱う。生体組織のような多孔質材料における浸透圧・流体流動などを解析。

  • 流体構造連成（FSI）解析
    流体と構造物の相互作用を解析可能。血流と血管壁などの連成挙動をモデル化。

  • 反応-拡散系解析
    化学反応と物質輸送を統合的に扱うことで、生体内の代謝反応や薬物動態モデルを構築可能。

バイオメカニクス特有の機能

  • 成長とリモデリング（Growth & Remodeling）
    組織の増殖、収縮、構造変化を時間的に追跡。血管の再構築や骨の成長などに応用可能。

  • 筋肉モデル
    筋収縮メカニズムを再現可能なアクティブマテリアルモデルを搭載。神経刺激への応答も設定可能。

  • 電気生理学的解析
    心筋電気伝導などを対象とした電気的信号の伝播解析が可能（FEBioと連携可能なFEBio Electric pluginが存在）。

  • 生体弁や血流モデル
    時間依存境界条件や動的荷重により、心臓弁や血流の流入・流出を再現。

モデル作成・可視化・連携機能

  • FEBio Studioとの統合
    GUIベースのモデリング、メッシュ作成、境界条件設定、結果可視化などを容易に行える。

  • メッシュ生成・編集機能
    六面体、四面体、シェル、ビームなど多様な要素を使用可能。メッシュの局所細分化も対応。

  • 境界条件・初期条件の柔軟な定義
    定数値、関数、表形式、時間依存の荷重や変位など、多様な条件を設定可能。

  • スクリプトやパラメトリック解析
    XMLベースの入力ファイルを用いて、パラメータスタディやバッチ処理の自動化が可能。

  • 外部データとの連携
    MRI・CTなどから取得した医用画像ベースのモデル作成が可能（Simplewareなどのツールと併用）。

  • 解析結果の可視化
    変形、応力、ひずみ、流速、濃度分布などの結果をカラーコンターやベクトル表示で視覚化可能。

FEBioの理論

FEBioでは構造解析と流体解析を行うことができる。 構造解析に関しては通常の支配方程式で表されるので省略する。 流体解析の支配方程式は以下のように表される。

連続方程式

$$ \begin{align} \dot{J} &= J \mathrm{div} \mathbf{v} \notag \end{align} $$

ここで、$J=\det\mathbf{F}$であり、ヤコビアン（現在時点と参照時点での体積比）を表す。 上付きのドットは時間微分を表す。 $\mathbf{F}$は変形勾配テンソル、$\mathbf{v}$は流速である。 この方程式はオイラー座標系で表現した連続方程式$\dot{\rho}+\rho \nabla\cdot\mathbf{v}=0$と等価である。

運動方程式

$$ \begin{align} \rho \mathbf{a} &= \mathrm{div} \mathbf{\sigma} + \rho \mathbf{b} \notag \end{align} $$

ここで、$\mathbf{a}=\frac{D\mathbf{v}}{Dt}$であり、加速度を表す。 $\rho$は密度、$\mathbf{\sigma}$はコーシー応力テンソル、$\mathbf{b}$は外力項である。 $\mathbf{\sigma}$は一般的に圧力の寄与(p)とせん断速度の寄与($\tau$)の２つの要因によって記述されることになる。

$$ \begin{align} \mathbf{\sigma} &= -p \mathbf{I} + \mathbf{\tau} \notag \end{align} $$

FEBioでは圧力$p$ではなく体積膨張率$e$を計算に使用する。 これらは次式で関係付けられる。

$$ \begin{align} p &= -K e \notag \end{align} $$

ここで$K$は体積膨張率(bulk modulus)である。

仮想仕事と弱形式

上記の支配方程式に含まれる未知変数は、体積膨張率$e$と流速$\mathbf{v}$である。 これらの未知変数を使いガレルキン有限要素法で離散化する。

$$ \begin{align} \delta W = \int_\Omega \delta \mathbf{v} \cdot (\mathrm{div} \mathbf{\sigma}+\rho(\mathbf{b}-\mathbf{a})) dv \notag \\ +\int_\Omega \delta J \left( \frac{\dot{J}}{J} -\mathrm{div} \mathbf{v}\right) dv \notag \end{align} $$

ここで、$\delta \mathbf{v}$は仮想速度、$\delta J$は仮想エネルギ密度に相当している。

発散定理を使い弱形式で書き直すと以下のように書き直す。

$$ \begin{align} & \delta W = \delta W_{ext} - \delta W_{int} = 0 \notag \\ & \delta W_{int} = \int_\Omega \mathbf{\tau} : \mathrm{grad} \delta \mathbf{v}dv +\int_\Omega \delta \mathbf{v} \cdot ( \mathrm{grad} p + \rho \mathbf{a} ) dv -\int_\Omega \left( \delta J \frac{\dot{J}}{J} + \mathrm{grad} \delta J \cdot \mathbf{v} \right) dv \notag \\ & \delta W_{ext} = \int_{\partial\Omega} \delta \mathbf{v} \cdot \mathbf{t}^{\mathbf{\tau}}da +\int_\Omega \delta \mathbf{v} \cdot \rho \mathbf{b} dv -\int_{\partial\Omega} \delta J v_n da \notag \end{align} $$

上記の方程式より自然境界条件は以下のように表されることがわかる。

  • 連続方程式：$v_n=0$ （境界における法線方向流速がゼロ)
  • 運動方程式：$\mathbf{t}^{\mathbf{\tau}}=0$ （接線方向の応力成分がゼロ）

流速の指定は、運動方程式の基本境界条件であるが間接的に連続方程式の自然境界条件にも 対応していることがわかる。

要素ライブラリ

以下のタイプのメッシュが有限要素メッシュとして利用可能である。 なお二次元モデルは使えないため、二次元モデルで計算する場合は厚さ方向に１要素をとった３次元メッシュを使用する。 要素の補間次数は１次と２次が利用可能である。

  • 6 面体要素 (brick element)
  • 5 面体要素 (wedge element)
  • 4 面体要素 (tetra element)

時間積分法

一般化α法を使用している。

FEBioの検証および妥当性

いくつかのVerification&Validationの計算結果を示す。 これらは、FEBioの信頼性を示す計算例になっている。 （なお計算の詳細については本記事の最後に示す参考URLを参照）

cavity cavity

step step

streching streching

billet billet

テスト計算例

例題を参考にして、2次元モデルを使い、超弾性体と液体の連成解析を行ってみた。

解析モデル

流入部を圧力指定（時間に関する波形パターンを指定）、流出部を圧力ゼロにし、流出入境界で変位をゼロとした。 流体が流入すると側面の弾性壁が変形するとした。 流体は水(密度$\rho=1000[kg/m^3]$、せん断粘度$\mu=0.001[Pa.s]$、体積弾性率$K=10^9[Pa]$) 、弾性壁は超弾性体材料 (Mooney-Rivlin材料：密度$\rho=1000[kg/m^3]$、体積弾性率$K=10^6[Pa]$、$c_1=1e4[Pa]$、$c_2=0[Pa]$)と設定した。

model

流入境界圧力の経時変化パターン（基準圧力=50[Pa]）

pattern

メッシュ分割

メッシュはgmshを使い作成した。 全解析領域を20節点6面体要素で分割し、abacus形式のフォーマット(inpファイル)で出力し、FEBioにインポートした。

mesh

解析結果

圧力コンターのアニメ

流速強度コンターのアニメ

渦度コンターのアニメ

まとめ

ここではオープンソースとして開発されている生体解析用プログラム FEBioを紹介した。 このソフトウエアの特徴として、以下のような利点が挙げられる。

  • 設定は統合環境内でGUIを介して行えるため、初心者にも取り付きやすい
  • 形状、メッシュなどを外部プログラムで作成してインポートすることができるため、使い慣れたツールを利用できる
  • 流体-構造連成解析が標準機能として組み込まれているため、流体解析と構造解析を個別に行い組み合わせる必要がない
  • ソースコードが公開されているためカスタマイズも可能

一方で、以下のような欠点もある。

  • 生体組織、生体材料の解析に特化しているため、複雑な熱輸送現象を扱うことができない
  • 乱流モデルを使えない
  • 現在のところ日本語で書かれた資料がほとんどない

実際に使ってみると統合環境GUI(FEBio studio)での設定は直感的で、また例題が多数公開されているのでそれほど手間を かけないで使うことができた。CAE解析（特に構造解析）のユーザは違和感なく使えるのではないかと感じる。 流体解析、構造解析を行っている担当者は必要に応じてFEBioを活用してみるとよいかもしれない。

参考資料1

参考資料2

弊社の解析事例

弊社の流体解析事例については、下記のリンクからご覧ください．

👉 解析事例のご紹介 FLUENT

👉 流体解析 実績及び事例