光を勉強するブログ

光と流体に関して。

ナビエ・ストークス方程式の導出(1)連続体について

導出の前に連続体という概念について記します。

 

1.連続体とは

 連続体、あるいは連続体近似とは、巨視的なスケールでの力学的挙動を表現するための数学的モデルのことであり、物理量を微視的スケールについて平均して得られる物理的性質を持つものを連続体といいます。ナビエ・ストークス方程式は流体の挙動を示す方程式ですが、例えば水の流れ、金属材料の変形といった挙動を考えるときにいちいち物質の原子・分子レベルまで考えていてはキリがありません。そこで巨視的に、つまり多数の原子・分子で構成される物質(系)で考えます。これは私達の感覚からしたら非常に小さいが、分子・原子の、いわゆる量子力学で扱うようなスケールからは遥かに大きいという一見都合の良い系になります。なので連続体では量子力学的な性質ではなくそれらを平均した結果あらわれる別の性質やその性質の時間変化を力学で記述することになり、これを巨視的に記述すると呼びます。また平均の結果表れる性質を巨視的性質、そこで用いられる物理量を巨視的物理量呼びます。そして、物質の持つ物理量が空間的に変化している場合でも、変化のスケールに対して微小な領域が定義でき、かつ巨視的な物理量が定義できるとき、その物体を連続体として扱えます。

 「微小領域」という言葉を用いましたが、これは十分に微小で空間的非一様性を無視できるが、それでも十分巨大で、その中に多数の原子分子が存在するようなものです。

 

2.連続体近似が妥当なための条件

 連続体近似が妥当であるためには、原子分子の力学が関与する時間・空間スケールに比べて十分に大きい微小領域が定義できる必要があり、その微小領域物体の振る舞いについて着目したい時間・空間スケールに比べ、十分に小さい局所的なものでなくてはなりません。

 

光拡散方程式の導出-(1)

 光拡散方程式は拡散光トモグラフィに用いられる重要な式であり、高散乱体内を一定距離進んだ位置において輸送方程式から近似的に導出されます。この近似を光拡散近似と呼ぶこともあるそうです。
 今回の記事から数回にかけて輸送方程式から光拡散方程式を導いてみようと思います。今回は目指すべきゴールの紹介と近似にあたって必要な準備を行います。

目標と方針

 目標は当然、光拡散方程式の導出です。そして導出を行う中での方針として、なるべく丁寧に式変形を追っていけたらと考えています。これは自分の理解を深め、今後より詳細な議論へ発展することができたらという思いがあることと、式変形を飛ばした導出はネット上にいくらでも転がっている、、、と思ったためです。

 では初めにスタート地点とゴール地点を紹介します。これは、あれこれ式変形して最後こうなりました!というものより、最後こうしたいからあれこれ変形していくよ!という方が自分好みだからです。

 前置きが長くなりましたが、ではまずスタート地点である輸送方程式を以下に示します。
\begin{align}
\displaystyle\frac{1}{c} \frac{\partial}{\partial t} I(\boldsymbol{r}, \boldsymbol{s}, t)+ \boldsymbol{s} \cdot \nabla I(\boldsymbol{r},\boldsymbol{s}, t)+(\mu_a +\mu_s)I(\boldsymbol{r}, \boldsymbol{s}, t) =\mu_s \int_{\boldsymbol{s}^2} p(\boldsymbol{s}',\boldsymbol{s})I(\boldsymbol{r},\boldsymbol{s'},t)d\boldsymbol{s'}\tag{1}
\end{align}

 この式の説明と導出は以前の記事を拝見していただけたらと思います。

 

iwalabo.hatenablog.com

 

 そして次にゴール地点です。(1)式をあれこれ変形して以下の式に近似していきたいと思います。
\begin{align}
\displaystyle \frac{1}{c} \frac{\partial}{\partial t}\Phi(\boldsymbol{r},t)=\nabla(D\nabla\Phi(\boldsymbol{r},t))-\mu_a\Phi(\boldsymbol{r},t) \tag{2}
\end{align}

 この式の詳細な説明はおいおい行いますが、ここで重要なことは光強度I(\boldsymbol{r}, \boldsymbol{s}, t)という位置、方向、時間の6つの変数を持った関数が、方向成分を失い\Phi(\boldsymbol{r}, t)という4の変数を持った関数になっているという点です。この関数は積分強度と呼ばれ、この関数を用いて近似を行っていきます。
 変数が2つも減っているので(1)式に比べてなんか解きやすそうですよね。

 

Beer Lambert則の歴史と導出

  近赤外分光法における重要な式、Beer Lambert則について書こうと思います。

歴史

  Lambert Beer則は名前の通り、Beerという人とLambertという人の功績からなっています。簡単に書くと、Lambertという人が、光強度が光源からの距離に対して指数関数的に減少することを発見し、Beerという人が媒質の濃度についても指数関数的に減少することを発見しました。これをまとめて、光が距離と濃度に指数関数的に減衰するという法則が定式化されました。
 また、Bougerという人がLambert以前に距離から指数関数的に減衰することを実験から示唆していたことからBouger Lambert Beer則とも呼ばれます。

 このことをまとめると、、、

f:id:Iwalabo:20180831113455p:plain

このようになります。始まりは3世紀も前なんですね。

ちなみにこのBougerは天文学に知見を持った方だったそうです。

 

導出

   ある媒体に光強度 I(\boldsymbol{r}, \boldsymbol{s}, t)で照射された光は、微小距離進む中で媒質による吸収の影響を受け、dIだけ変化します。
 また、I_{吸収}\Delta lの関係は以下のような式で記述されることが知られています。

\begin{align}
\displaystyle \Delta I_{吸収}=-\mu_a I \Delta l \tag{1}
\end{align}

 これは、吸収による光強度の変化量が変化前の光強度に比例し、その比例定数が吸収係数\mu_a[mm^{-1}]であることを示しています。そして符号が負なのは吸収によって光強度が減少していることを示しています。
 また、ここでΔを無視できるくらいの微小量とすることで、(1)式は

\begin{align}
\displaystyle dI &= -\mu_a I dl \tag{2} \\ \\
\displaystyle \frac{dI}{dl} &= -\mu_a I \tag{3}
\end{align}

となります。

 (3)式のようにすることで簡単な微分方程式となりますのでこれを解いていきます。
\begin{align}
\displaystyle \frac{1}{I}dI &= -\mu_a dl \\ \\
\displaystyle \int \frac{1}{I}dI &= -\mu_a \int dl \\ \\
\displaystyle logI+C_1 &= -\mu_a l + C_2 \\ \\
\displaystyle I &= \exp(-\mu_a l+C_3) \\ \\
\displaystyle I &= C_4 \exp(-\mu_a l) \\ \\
\end{align}
  これで一般解が求まりました。ここでC_1,C_2,C_3,C_4は0以外の実数です。
次にl=0においてI = I_0という境界条件のもと特殊解を求めます。
\begin{align}
\displaystyle I_0 &= C \exp(-\mu_a \times 0) \\ \\
\displaystyle C &= I_0 \\ \\
\end{align}
よって特殊解は
\begin{align}
\displaystyle I &= I_0 \exp(-\mu_a l) \tag{4}
\end{align}
となります。
  ここで吸収係数\mu_a = \varepsilon Cという関係を用います。ここで\varepsilonはモル吸光係数という物質固有の値で、Cは媒質の濃度になります。すると(4)式は
\begin{align}
\displaystyle I &= I_0 \exp(-\varepsilon C l) \tag{5}
\end{align}
となり、光強度が光源からの距離と媒質の濃度に対して指数関数的に減衰するという式が導出されました。これがLambert Beer則になります。

 図にすると以下のようなイメージです。   

            f:id:Iwalabo:20180831123457p:plain

光の伝搬現象における輸送方程式の概要と導出

概要と方針

 輸送方程式(Radiative Transfer Equation : RTE)とはボルツマン方程式の一種である偏微分方程式です。光計測における基礎方程式として知られ、光CTや拡散光トモグラフィに利用されていものです。

   さて、この記事の方針ですが、題目の通り輸送方程式を光輸送の観点から解析学に知見の浅い方でも分かるよう導出していきます。というのも、以前、ネットで検索したら中性子関連のものだらけで光に関するものが(あまり)見つからなく、あったとしても難解なものだらけでした。初学者であり、解析学に知見の浅い自分でも分かるような、直観でイメージできる記事があったら…と考え、ないなら自分で書いてしまおうと思ったことが動機です。

 

導出

 輸送方程式はエネルギー保存則を考えることで導くことが叶います。これは光の伝搬現象を導く際にも同様です。具体的には光強度の空間的・時間的な変化は光の吸収と散乱によっておこる。ということを記述します。式で表すと以下のようになります。

\begin{align}
\displaystyle I(\boldsymbol{r}+\Delta l\boldsymbol{s}, \boldsymbol{s}, t+\Delta t) = I(\boldsymbol{r}, \boldsymbol{s}, t) + \Delta I_{吸収} + \Delta I_{散乱}  \tag{1}
\end{align}

  {I(\boldsymbol{r}, \boldsymbol{s}, t)}はある位置\boldsymbol{r} \in \boldsymbol{R}^3、方向 \boldsymbol{s}\in \boldsymbol{S}^2、時刻tにおける光強度を表しており、 \Delta I_{吸収} + \Delta I_{散乱}はそれぞれ吸収と散乱による光強度の変化を表します。この式の左辺では\boldsymbol{s}方向に微小距離\Delta l移動し、時刻\Delta t経過したときの光強度の微小変化を示し、右辺で、その変化量は吸収と散乱によると表しています。まず、この左辺についてテイラー展開をもちいて一次近似を施します。すると

\begin{align}
\displaystyle L.H.S &=  I(\boldsymbol{r}+\Delta l\boldsymbol{s}, \boldsymbol{s}, t+\Delta t) \\ \displaystyle &= \sum_{k=0}^1 \frac{1}{k!} (\Delta l \boldsymbol{s} \cdot \nabla + \Delta t \frac{\partial}{\partial t} )^k I(\boldsymbol{r}, \boldsymbol{s}, t) \\ \displaystyle  &= I(\boldsymbol{r}, \boldsymbol{s}, t)+\Delta l \boldsymbol{s} \cdot \nabla  I(\boldsymbol{r}, \boldsymbol{s}, t) + \Delta t \frac{\partial}{\partial t}I(\boldsymbol{r}, \boldsymbol{s}, t) \tag{2}
\end{align}

そして(1)式右辺の I(\boldsymbol{r}, \boldsymbol{s}, t)を導出された(2)式左辺に移項することによって(1)式は次のようになります。

\begin{align}
\displaystyle \Delta l \boldsymbol{s} \cdot \nabla  I(\boldsymbol{r}, \boldsymbol{s}, t) + \Delta t \frac{\partial}{\partial t}I(\boldsymbol{r}, \boldsymbol{s}, t) = \Delta I_{吸収} + \Delta I_{散乱}  \tag{3}
\end{align}

 さて、次に光の吸収による光強度の変化I_{吸収}をみます。

               f:id:Iwalabo:20180830112043p:plain

上図のように、ある媒体に光強度 I(\boldsymbol{r}, \boldsymbol{s}, t)で照射された光は、微小距離進む中で媒質による吸収の影響を受け、dIだけ変化します。
 また、I_{吸収}\Delta lの関係は以下のような式で記述されることが知られています。

\begin{align}
\displaystyle \Delta I_{吸収}=-\mu_a I(\boldsymbol{r}, \boldsymbol{s}, t) \Delta l \tag{4}
\end{align}

 これは、吸収による光強度の変化量が変化前の光強度に比例し、その比例定数が吸収係数\mu_a[mm^{-1}]であることを示しています。そして符号が負なのは吸収によって光強度が減少していることを示しています。

…と書きましたがこの(4)は正確ではありません。

 正確にはI_{吸収}dIとおき、\Delta ldlとおきかえることが叶うときに成り立ちます。また、そうすることでこの(4)は単純な微分方程式となり、変数分離型として解くことで光トポグラフィにおける重要な式、Beer Lambert則を得ることができます。

 

iwalabo.hatenablog.com

 

 

 次に散乱による光強度の変化I_{散乱}をみます。散乱の影響は吸収によるものと比べて少し複雑になります。また、輸送方程式において散乱による項は2つ現れます。散乱による光強度の増加項と散乱による光強度の減少光の2つです。

              f:id:Iwalabo:20180830112114p:plain

  Fig.2は散乱によって生ずる光強度の減少を表しています。\boldsymbol{s}方向からの光の一部が媒質内で粒子に衝突し、\boldsymbol{s}’方向へと経路を変えることで本来の\boldsymbol{s}方向の光が減少します。この現象は以下のように記述されます。

\begin{align}
\displaystyle \Delta I_{散乱-}=-\mu_s \int_{\boldsymbol{s}^2} p(\boldsymbol{s},\boldsymbol{s}')I(\boldsymbol{r},\boldsymbol{s},t)d\boldsymbol{s}' \Delta l \tag{5}
\end{align}

 ここで、符号の負は先ほどと同様に光強度の減少を示し、\mu_sは散乱係数[mm^{-1}]を示します。そして、p(\boldsymbol{s},\boldsymbol{s}')は方向\boldsymbol{s}から方向\boldsymbol{s}'への散乱確率を示す位相関数になります。この位相関数と光強度の積をとり、\boldsymbol{s}'に関して全立体角で積分することで\boldsymbol{s}方向の光強度が別方向にどれだけ散乱されるかわかります。また、位相関数p(\boldsymbol{s},\boldsymbol{s}')は全立体角で積分すると1になるよう定義されており、これを用いてこの積分式を整理すると、

\begin{align}
\displaystyle \Delta I_{散乱-} &= -\mu_sI(\boldsymbol{r},\boldsymbol{s},t) \int_{\boldsymbol{s}^2} p(\boldsymbol{s},\boldsymbol{s}')d\boldsymbol{s}' \Delta l \\ \displaystyle &= -\mu_s I(\boldsymbol{r},\boldsymbol{s},t)\Delta l \tag{6}
\end{align}

 次に散乱による光強度の増加をみます。 

                  f:id:Iwalabo:20180830114246p:plain

  Fig.3は散乱によって生ずる光強度の増加を表しています。減少によるものと同様に、\boldsymbol{s}'方向からの光の一部が媒質内で粒子に衝突し、\boldsymbol{s}方向へと経路を変えることで本来の\boldsymbol{s}方向の光が増加します。この現象は以下のようになります。

\begin{align}
\displaystyle \Delta I_{散乱-}=\mu_s \int_{\boldsymbol{s}^2} p(\boldsymbol{s}',\boldsymbol{s})I(\boldsymbol{r},\boldsymbol{s'},t)d\boldsymbol{s'} \Delta l \tag{7}
\end{align}
  先ほどと異なるのは位相関数の変数の順序と、光強度Iの方向変数が\boldsymbol{s}'となっており、増加現象を示しているため符号が正となっている点です。

 最後に、これまで導出した吸収・散乱項(4)、(6)、(7)式を(3)式に代入すると、
\begin{align}
\displaystyle \Delta l \boldsymbol{s} \cdot \nabla I(\boldsymbol{r},\boldsymbol{s}, t) + \Delta t \frac{\partial}{\partial t} I(\boldsymbol{r}, \boldsymbol{s}, t) =  -\mu_a I(\boldsymbol{r}, \boldsymbol{s}, t) \Delta l -\mu_s I(\boldsymbol{r},\boldsymbol{s},t)\Delta l +\mu_s \int_{\boldsymbol{s}^2} p(\boldsymbol{s}',\boldsymbol{s})I(\boldsymbol{r},\boldsymbol{s'},t)d\boldsymbol{s'} \Delta l \tag{8}
\end{align}
となります。そして全項を\Delta lで割り、媒質中の光速度c=\frac{\Delta l}{\Delta t}の関係式を用いると、
\begin{align}
\displaystyle \boldsymbol{s} \cdot \nabla I(\boldsymbol{r},\boldsymbol{s}, t) + \frac{1}{c} \frac{\partial}{\partial t} I(\boldsymbol{r}, \boldsymbol{s}, t) =  -\mu_a I(\boldsymbol{r}, \boldsymbol{s}, t)-\mu_s I(\boldsymbol{r},\boldsymbol{s},t)+\mu_s \int_{\boldsymbol{s}^2} p(\boldsymbol{s}',\boldsymbol{s})I(\boldsymbol{r},\boldsymbol{s'},t)d\boldsymbol{s'}\tag{9}
\end{align}
となり、整理することで輸送方程式を得ることができます。
\begin{align}
\displaystyle\frac{1}{c} \frac{\partial}{\partial t} I(\boldsymbol{r}, \boldsymbol{s}, t)+ \boldsymbol{s} \cdot \nabla I(\boldsymbol{r},\boldsymbol{s}, t)+(\mu_a +\mu_s)I(\boldsymbol{r}, \boldsymbol{s}, t) =\mu_s \int_{\boldsymbol{s}^2} p(\boldsymbol{s}',\boldsymbol{s})I(\boldsymbol{r},\boldsymbol{s'},t)d\boldsymbol{s'}\tag{10}
\end{align}

 

  さて、輸送方程式を導くことが出来ました。なるべく直感的にわかりやすく、かつ曖昧な議論を避けたつもりですが、良く分からない点や間違いがあったら申し訳ありません。

  今後は輸送方程式を用いて光拡散近似を行っていこうと考えております。興味があったらご一読お願いします。

 

iwalabo.hatenablog.com