基礎方程式の空間差分近似(空間についての離散化) Taylor展開の便利さと言ったら・・・

非圧縮性流体

方針

一旦、熱はややこしくなるので考えずに、どうやったらNavier-Stocks方程式が解けるのか考えていきましょう。
Navier-Stocks方程式は以下のように表せるのでした。

$$
\begin{align}
& \frac{ \partial u_i }{ \partial x_i } = 0 \tag{1} \label{eq:scale-factor-1} \\\\
& \frac{ \partial u_i }{ \partial t} + \frac{ \partial u_i u_j }{ \partial x_j} = -\frac{1}{\rho} \frac{\partial P}{\partial x_i} + \nu \frac{\partial^2 u_i}{\partial x_j^2} \tag{2} \label{eq:scale-factor-2} \qquad
\end{align}
$$

ベクトル式

$$
\begin{align}
&\nabla \cdot \boldsymbol V = 0 \tag{1′} \label{eq:scale-factor-101} \\\\
&\frac{ \partial \boldsymbol V}{ \partial t} + (\boldsymbol V \cdot \nabla) \boldsymbol V = -\frac{1}{\rho} \nabla p + \nu \nabla^2 \boldsymbol V \label{eq:scale-factor-102} \qquad \\\\
\end{align}
$$

方針としては2つです。
① \(\int_t \int_V \)のように積分すれば解けるはず?
② \eqref{eq:scale-factor-1}と\eqref{eq:scale-factor-2}を連立して解く。

①に関しては、単純に考えて、これは時間と空間に関する偏微分方程式なので、時間と空間で積分してやれば解けるはずです。しかし、まあそれができたら誰も苦労しません。Navier-Stocks方程式は解析的に解けるどころか、滑らかな解の存在さえ、証明されていません。方程式の 3次元空間と(1次元の)時間の中で初期速度を与えると、ナビエ–ストークス方程式の解となる速度ベクトル場と圧力のスカラー場が存在して、双方とも滑らかで大域的に定義されるだろうか、という問題がミレニアム懸賞問題となっています。ちゃんと滑らかな解があるかすら分かっていません。(離散的な解は数値解析的に解けるので存在はしているが、滑らかな解があるかはわからない。また、2次元の限定的な場面では解析解が得られる。)ミレニアム懸賞問題とは、数学的に難問かつ重要な問題に対し、アメリカのクレイ数学研究所が100万ドルの懸賞金を掛けている問題です。勇気を持って人生を捨てれる人はぜひこの未解決問題に取り組んでみて欲しいと思います。

移流項が難しいです。移流項が。移流項は、\(u_i u_j\)なので、言ってみれば2乗の項です。2乘の項はいわゆる非線形問題というやつです。例えば、\(\frac{\partial u^2}{\partial x}\)のuに、sin成分が入ってくるとします。その時、

$$
\begin{align}
\frac{ \partial u^2 }{ \partial x } &= \frac{ \partial \sin^2 x }{ \partial x } \\ &= \frac{ \partial \sin x}{\partial x} \sin x + \sin x \frac{\partial \sin x}{ \partial x }  \\ &= \cos x \sin x + \sin x \cos x \\ &= \sin 2x \qquad \tag{3} \label{eq:scale-factor-3}
\end{align}
$$

これが何を意味しているかというと、微分されて出ていくものの周波数が2倍になっています。どんどん細かくなっていきます。そして、フーリエ級数展開をより、だいたいの関数は、sin成分に分解できるので、周期成分が時間を経るごとに細かくなっていきます。流体的なイメージとしては、渦が分解されていきます。大きな渦が分解されていき、より細かな渦になっていく。これが乱流であり、やっかいなものなのです。解く難易度が一気に上がります。これが根本的な難しさです。

②のように\eqref{eq:scale-factor-1}の連続の式をどうやって反映させるかも難しいところです。今は未知数は\(u_1, u_2, u_3, P\)の4つで、\eqref{eq:scale-factor-2}のNavier-Stocks方程式は3つの式だから、連続の式を使ってやらないと解けません。どうやって連続の式を使うか、難しいところです。これについては数値解析的にはうまく使用する方法があり、おそらく感動することになると思います。


Taylor展開

ではどうすればいいのでしょうか?微分で解けないなら差分にするしかありません。みんな大好きTaylor展開です。そして離散化します。
Taylor展開の前に、微分の定義から攻めましょう。\(\Delta x\)が十分に小さい時、

$$
\begin{align}
\left.\frac{\partial f}{\partial x}\right\vert_x &= \lim_{ \Delta x \to 0} \frac{ f(x + \Delta x) – f(x) }{ \Delta x } \\
&\simeq \frac{ f(x + \Delta x) – f(x) }{ \Delta x } \\\\
f(x &+ \Delta x) \simeq  f(x) +  \Delta x \left.\frac{\partial f}{\partial x}\right\vert_x \tag{4} \label{eq:scale-factor-4}
\end{align}
$$

のように近似することができます。差分近似と言います。そして、これはTaylor展開の第1次近似とも言えます。
Taylor展開と言うのはなんでしたでしょうか。

$$
\begin{align}
f(x + \Delta x) &= f(x)+\frac{\Delta x}{1!} \left.\frac{\partial f}{\partial x}\right\vert_x +\frac{\Delta x^2}{2!} \left.\frac{\partial^2 f}{\partial x^2}\right\vert_x +\frac{\Delta x^3}{3!} \left.\frac{\partial^3 f}{\partial x^3}\right\vert_x + \cdots  \tag{5} \label{eq:scale-factor-5}
\end{align}
$$

これです。簡単に証明と解説もしておきます。勉強になるので知っておいて損はないでしょう。要は、xから少し離れたf(x)を、f(x)の値とf(x)の微分値を使って近似的に求めてやろうと言う方法です。Taylor展開の2階微分の項以降を捨て去った形が、微分の定義式からの差分近似です。
そして、もう少し数値解析的に考えます。x軸のある範囲を均等にN個の領域で分割したと考えて、分割した点にそれぞれ1からn+1まで番号を振っていきます。この時、i点周りのTaylor展開で、i番目の点の値\(f_i\)とその微分値から、\(f_{i+1}\)の点の値を求めることができます。つまり、

$$
\begin{align}
f_{i+1} &= f_i+\frac{\Delta x}{1!} \left.\frac{\partial f}{\partial x}\right\vert_i +\frac{\Delta x^2}{2!} \left.\frac{\partial^2 f}{\partial x^2}\right\vert_i +\frac{\Delta x^3}{3!} \left.\frac{\partial^3 f}{\partial x^3}\right\vert_i + \cdots  \tag{6} \label{eq:scale-factor-6}
\end{align}
$$

同様にして、

$$
\begin{align}
f_{i-1} &= f_i-\frac{\Delta x}{1!} \left.\frac{\partial f}{\partial x}\right\vert_i +\frac{\Delta x^2}{2!} \left.\frac{\partial^2 f}{\partial x^2}\right\vert_i -\frac{\Delta x^3}{3!} \left.\frac{\partial^3 f}{\partial x^3}\right\vert_i + \cdots  \tag{7} \label{eq:scale-factor-7}
\end{align}
$$

ともできます。ここで、Navier-Stocksは微分方程式なので、微分値が欲しい。つまり、\(f_i\)や\(f_{i+1}\)の値から、逆に微分値を求めていくと言うことをします。\eqref{eq:scale-factor-6}から、

$$
\begin{align}
\left.\frac{\partial f}{\partial x}\right\vert_x &= \frac{f_{i+1} – f_{i}}{\Delta x}-\frac{\Delta x}{2} \left.\frac{\partial^2 f}{\partial x^2}\right\vert_x + \cdots  \tag{8} \label{eq:scale-factor-8}
\end{align}
$$

そして、2項目以降を打ち切ることになります。この時の誤差の第1項目が\(\Delta x\)に比例しているので、第1次近似、1次精度差分近似と言われます。

ここで、\(f_{i+1} – f_{i-1}\)を考えてみる。すると、

$$
\begin{align}
&f_{i+1} – f_{i-1} = \frac{2 \Delta x}{1!} \left.\frac{\partial f}{\partial x}\right\vert_i + \frac{2 \Delta x^3}{3!} \left.\frac{\partial^3 f}{\partial x^3}\right\vert_i + \cdots  \tag{9} \label{eq:scale-factor-9} \\\\ &\left.\frac{\partial f}{\partial x}\right\vert_i =\frac{f_{i+1} – f_{i-1}}{2\Delta x} – \frac{\Delta x^2}{6} \left.\frac{\partial^3 f}{\partial x^3}\right\vert_i + \cdots  \tag{10} \label{eq:scale-factor-10}
\end{align}
$$

となります。これは、誤差の第1項目が\(\Delta x^2\)に比例しています。つまり、先ほどより打ち切り誤差が少ないと言えます。これは、2次精度差分近似になります。すげー、うひゃひゃひゃ。このように微分値を近似的に求めてやり、それを微分方程式にぶっこんで差分方程式にしてしまえば方程式が解けてしまいます。全ての点で式を一つづつ立てて、連立方程式にして解くと言う感じです。両端の値はもちろん与えてやらなければならなく、これを境界条件と言ったりします。かっちょいい。

同様に、3次精度、4字精度差分近似と言うものを考えることができる。長くなるけどやってみよす。


3次精度・4次精度差分近似

同じように、i点周りのTaylor展開を考えるのだが、2つしか式がないと、2次精度が限界なので、2個分離れた場所のTaylor展開を考えていきます。\(\left.\frac{\partial }{\partial }\right\vert_i\)を書くのが少し面倒なので省略します。\(\frac{\partial f}{\partial x}\)は、i点の微分値だと思ってください。

$$
\begin{align}
&f_{i+2} = f_i+\frac{2\Delta x}{1!} \frac{\partial f}{\partial x} +\frac{(2\Delta x)^2}{2!} \frac{\partial^2 f}{\partial x^2} +\frac{(2\Delta x)^3}{3!} \frac{\partial^3 f}{\partial x^3} +\frac{(2\Delta x)^4}{4!} \frac{\partial^4 f}{\partial x^4}+ \cdots  \tag{11} \label{eq:scale-factor-11} \\\\
&f_{i+1} = f_i+\frac{\Delta x}{1!} \frac{\partial f}{\partial x} +\frac{\Delta x^2}{2!} \frac{\partial^2 f}{\partial x^2} +\frac{\Delta x^3}{3!} \frac{\partial^3 f}{\partial x^3} +\frac{\Delta x^4}{4!} \frac{\partial^4 f}{\partial x^4}+ \cdots  \tag{12} \label{eq:scale-factor-12} \\\\
&f_{i-1} = f_i-\frac{\Delta x}{1!} \frac{\partial f}{\partial x} +\frac{\Delta x^2}{2!} \frac{\partial^2 f}{\partial x^2} -\frac{\Delta x^3}{3!} \frac{\partial^3 f}{\partial x^3} +\frac{\Delta x^4}{4!} \frac{\partial^4 f}{\partial x^4}- \cdots  \tag{13} \label{eq:scale-factor-13} \\\\
&f_{i-2} = f_i-\frac{2\Delta x}{1!} \frac{\partial f}{\partial x} +\frac{(2\Delta x)^2}{2!} \frac{\partial^2 f}{\partial x^2} -\frac{(2\Delta x)^3}{3!} \frac{\partial^3 f}{\partial x^3} +\frac{(2\Delta x)^4}{4!} \frac{\partial^4 f}{\partial x^4}- \cdots  \tag{14} \label{eq:scale-factor-14}
\end{align}
$$

ここで、3次精度差分は、前3つを、4次精度差分は全て使って、微分項を消していきます。条件を求めるため、係数を表にしてみます。

 \(\frac{\partial f}{\partial x}\)\(\frac{\partial^2 f}{\partial x^2}\)\(\frac{\partial^3 f}{\partial x^3}\)\(\frac{\partial f}{\partial x}\)\(\frac{\partial^2 f}{\partial x^2}\)\(\frac{\partial^3 f}{\partial x^4}\)\(\frac{\partial^4 f}{\partial x^4}\)
\(f_{i+2}\)2\(\frac{4}{2}\)\(\frac{8}{6}\)2\(\frac{4}{2}\)\(\frac{8}{6}\)\(\frac{16}{24}\)
\(f_{i+1}\)1\(\frac{1}{2}\)\(\frac{1}{6}\)1\(\frac{1}{2}\)\(\frac{1}{6}\)\(\frac{1}{24}\)
\(f_{i-1}\)-1\(\frac{1}{2}\)\(-\frac{1}{6}\)-1\(\frac{1}{2}\)\(-\frac{1}{6}\)\(\frac{1}{24}\)
\(f_{i-2}\)-2\(\frac{4}{2}\)\(-\frac{8}{6}\)\(\frac{16}{24}\)
1001000

それぞれの式をa、b、c、d倍して足し合わせるとして、表の縦の係数を足して、\(\frac{\partial f}{\partial x}\)が1になって、他は0になるように足してやれば良いです。即ち、

\( \begin{eqnarray} \left\{ \begin{array}{l} 2a + b -c = 1 \\ \frac{4}{2}a + \frac{1}{2}b + \frac{1}{2}c = 0 \\ \frac{8}{6}a + \frac{1}{6}b – \frac{1}{6}c = 0 \end{array} \qquad \Leftrightarrow \qquad \right.\left\{ \begin{array}{l} a =-\frac{1}{6} \\ b = 1 \\ c = -\frac{1}{3} \end{array} \right. \end{eqnarray} \)

\( \begin{eqnarray} \left\{ \begin{array}{l} 2a + b -c -2d = 1 \\ \frac{4}{2}a + \frac{1}{2}b + \frac{1}{2}c + \frac{4}{2}d = 0 \\ \frac{8}{6}a + \frac{1}{6}b – \frac{1}{6}c – \frac{8}{6}d= 0 \\ \frac{16}{24}a + \frac{1}{24}b + \frac{1}{24}c + \frac{16}{24}d= 0 \end{array} \qquad \Leftrightarrow \qquad \right.\left\{ \begin{array}{l} a =-\frac{1}{8} \\ b = 1 \\ c = -1 \\ d=\frac{1}{8} \end{array} \right. \end{eqnarray} \)

よって、3次精度差分は、

$$
\begin{align}
&- \frac{1}{6} f_{i+2} + f_{i+1} – \frac{1}{3} f_{i-1} = \frac{1}{2} f_i + \Delta x \frac{\partial f}{\partial x} + O(\Delta x^4) \\\\
&\Leftrightarrow \frac{\partial f}{\partial x} = \frac{- f_{i+2} + 6f_{i+1} – 3f_i – 2f_{i-1}}{6 \Delta x} + O(\Delta x^3) \tag{15} \label{eq:scale-factor-15}
\end{align}
$$

同じように、4次精度差分は、

$$
\begin{align}
&- \frac{1}{8} f_{i+2} + f_{i+1} – f_{i-1} + \frac{1}{8} f_{i+2} = \frac{8}{12} \Delta x \frac{\partial f}{\partial x}+ O(\Delta x^5) \\\\
&\Leftrightarrow \frac{\partial f}{\partial x} = \frac{- f_{i+2} + 8f_{i+1} – 8f_{i-1} + f_{i-2}}{12 \Delta x} + O(\Delta x^4)\tag{16} \label{eq:scale-factor-16}
\end{align}
$$

よっしゃーこれで精度が高い計算ができるぜーと思いきや、そうも行きません。見ての通り、2個以上離れた点の情報を使ってしまっています。こうなると、端の隣の格子点になると適用できません。何故なら、境界条件よりも外の点の値など定義できないからです。安易に境界条件同じ値を使うなどと考えては行けません。境界条件はシビアであり、固定条件以外にも色々なものがあります。よって、これらは内側の部分でしか使えないのです。これはあまり望ましくありません。下の図を見てください。

まずは、Taylor展開で差分近似し、で次の点の値を求めるとしましょう。図はわかりやすいように1次精度差分近似としました。本来の微分値が求まり、その点の値と微分値からTaylor展開で次の点の値を求めます。Taylor展開と言うのは、1次近似の場合その点から( \(傾き\times \Delta x\))で次の点を求めるのでした。図1 を見てもらうとわかるが、変化率が大きいほど近似する時に誤差が生まれやすいです。また、今回は微分値を、差分近似して求めるのでした。この話は図2のようになります。本来の微分値と差分近似して求めた微分値は、こちらも変化率が大きいほど誤差が生まれてしまいます。普通に物理現象を考えると最も変化率が大きいのは境界付近です。境界から物理量が伝わって変化していくのですから。だいたい、微分方程式系の数値解析は初期値・境界値問題なのです。うわ〜。よって、最も精度が欲しい境界付近に適用できないなら高次精度ってあまり意味がないんじゃ・・・と思ってしまいます。

それに、参照点も増えて計算量が単純に増えます。なら、メッシュを細かくすりゃいいじゃん。メッシュを細かくすれば問答無用で\(\Delta x\)が小さくなって誤差が大きく減る。境界条件付近も問題ないし。計算量増えるけど、どうせ増えるならそっちの方が良くね?内側だけに適用するのもめんどくさいし・・・と言うことで、精度を高める時、真っ先にメッシュを細かくするという方法が適用がされます。高次精度にするのはその後です。

これらの差分近似を、流速に対しては、1次精度風上差分、2次精度中心差分、3次精度風上差分、4次精度風上差分と言います。奇数精度のものは”風上差分”です。流速というのは、流速が速い方から遅い方へ伝わっていきます。よって、遅い方の情報から速い方の値を求める解くことはできません。めちゃくちゃな値しか出ないと思います。なので、奇数精度の場合は、必ず風上から解いていくので、風上差分と言う名前が付いています。

本来の意味からすると、奇数精度の風上差分が良いと思われるのだが、適用する項によっても違います。移流項に対しては、運ばれていくような物理現象なので、風上差分がいいのですが、拡散項、圧力項に関しては中心差分の方が良さそうな気がします。これは物理現象からも、高い方から低い方と言ったものではなく双方向で、周りに馴染んでいくよう物理量だからです。そして、一つの項だけ、奇数精度でやるのは、ちょっと解きにくかったりして、実際には移流項も中心差分でやってしまうことも多かったりします。

 


<

p style=”text-align: center;”>
[
前の記事へ]  [非圧縮性流体の目次へ]  [次の記事へ]