数学とか語学とか楽しいよね

フランス語、ドイツ語、ロシア語、オランダ語、英語、スペイン語、ラテン語とか数学とか数値計算(有限要素法、有限体積法、差分法、格子ボルツマン法、数理最適化、C++コード付き)とか勉強したことをまとめます。右のカテゴリーから興味のある記事を探してください。最近はクラシックの名演も紹介しています。noteにも書いています。https://note.mu/kumag_min

【有限要素法】1次元非定常移流拡散方程式の解析解を導く

今回は一次元非定常移流拡散方程式の解析解を導きたいと思います(実は途中まで)。有限要素法のプログラムがきちんと動くかどうか確かめるために使いたいのです。これ意外と本に載っていないのではないでしょうか。

一次元非定常移流拡散方程式とは

 \displaystyle \frac{\partial u}{\partial t}+V\frac{\partial u}{\partial x}-D\frac{\partial^2 u}{\partial x^2}=0 in  (-\infty, \infty)

のような無限領域で定義された偏微分方程式のこととし、初期条件はある関数  f(x) で与えられるとします。ここで、 u(t,x) は未知数で時間と空間の二変数関数、 V は流速で定数、 D は拡散係数で定数です。熱またはある溶質が時間が経過するにつれてどのような分布をとるか(移流拡散していくか)求めるのがこの問題です。無限領域なので境界条件は必要なく、初期条件だけです。

まず、以下のように座標を  (t,x) から  (\tau,y) へと変換します。


 
\displaystyle 
\begin{eqnarray} 

\tau &=& t \\
y &=& x-Vt \\

\end{eqnarray}

それぞれの偏微分を計算しておきましょう。

 \displaystyle \frac{\partial \tau}{\partial x}=0, \frac{\partial \tau}{\partial t} =1, \frac{\partial y}{\partial t}=-V, \frac{\partial y}{\partial x}=1

このとき空間と時間の一階微分はそれぞれ



\displaystyle 
\begin{eqnarray} 

\frac{\partial u}{\partial x} &=& \frac{\partial u}{\partial \tau}\frac{\partial \tau}{\partial x} + \frac{\partial u}{\partial y}\frac{\partial y}{\partial x} \\
 &=& 0 +  \frac{\partial u}{\partial y}\\

\frac{\partial u}{\partial t} &=& \frac{\partial u}{\partial \tau}\frac{\partial \tau}{\partial t} + \frac{\partial u}{\partial y}\frac{\partial t}{\partial x} \\
 &=& \frac{\partial u}{\partial \tau} -  V\frac{\partial u}{\partial y}\\

\end{eqnarray}

と計算できます。これらをもとの式に代入すると



\displaystyle 
\begin{eqnarray} 

\frac{\partial u}{\partial t}+V\frac{\partial u}{\partial x}-D\frac{\partial^2 u}{\partial x^2} &=& \frac{\partial u}{\partial \tau} -  V\frac{\partial u}{\partial y} +V\frac{\partial u}{\partial y} -D\frac{\partial^2 u}{\partial y^2}\\
 &=& \frac{\partial u}{\partial \tau} -D\frac{\partial^2 u}{\partial y^2} =0

\end{eqnarray}

となります。すなわち

 \displaystyle \frac{\partial u}{\partial \tau} = D\frac{\partial^2 u}{\partial y^2}

であり、最終的に拡散方程式に変換できるのです!初期条件は  f(y+Vt)=f(y) です。これはフーリエ変換を使って簡単に解くことができます。フーリエ変換 u(\tau,y) について解いて、座標の逆変換で  u(t,x-Vt) としてやればもとの非定常移流拡散方程式の解析解が得られます。

上記の解は

 \displaystyle u(\tau,y) = \frac{1}{\sqrt{4D\tau\pi}} \int_{-\infty}^{\infty} f(s)  e^{-\frac{(y-s)^2}{4D\tau}} ds

となります。なので、座標の逆変換  u(t,x-Vt) を使うと解の表示が

 \displaystyle u(t,x) = \frac{1}{\sqrt{4Dt\pi}} \int_{-\infty}^{\infty} f(s)  e^{-\frac{(x-Vt-s)^2}{4Dt}} ds

のように得られます。解析解と初めに言いましたが、この積分を解析的に求められる場合以外は解析解にはならないのでご注意下さい。ただ  f(x) をうまく選べば大丈夫です。