トップ > ライブラリ > 数値解析 > ラプラス変換 > ラプラス変換・例題５

数値解析：ラプラス変換・例題５

　現在工事中です。
　ラプラス変換を適用し、工学分野で現れる偏微分方程式を解くことができる。その具体例を 以下に紹介する。ここで示す解析の例題５では、数値解析・ラプラス変換にて解説したラプラス変換の概要を踏まえ、 解析解の式を具体的に導出している。
　本ページでは例題５の、直交座標系の半無限固体中の非定常伝熱問題を取り扱っているが、
　例題６では直交座標系の粒子充填層（拡散モデル）の滞留時間分布の問題を取り扱う。
　例題７では円筒座標系の非定常拡散問題を取り扱う。

先頭に戻る

【注】数式表示にはMathJaxを利用しています。IE8以下では表示が遅くなる可能性があります。FireFox などIE8以外のブラウザを利用下さい。

半無限固体中の非定常伝熱（直交座標系）

　半無限固体中の非定常伝熱方程式の解を求める例題５の解法を以下に示す。
　例題を再掲すると、次のとおり。

例第５：半無限固体中の非定常伝熱5)

　\(z=0\ \)から\(z=\infty\)まである半無限の固体が最初\(T=T_i\)に保たれ、時間\(t \gt 0\)で \(z=0\)の表面温度が\(T_w\)となりそのまま維持されるとき、固体中の温度分布の時間変化を求めよ。
　z方向の熱伝導の基礎式(5-1)式、および初期条件(5-2)式、境界条件(5-3)式は以下のとおりとする。熱伝導度\(\lambda\)は定数とす。 \[ \begin{align*} \frac{\partial T}{\partial t}=\lambda \frac{\partial^2 T}{\partial z^2} \tag{5-1} \end{align*} \] ここで、\(\lambda\)は固体の熱伝導度、\(T\)は温度、\(t\)は時間、\(z\)は表面からの距離とする。 \[ \begin{align*} {\rm I.C.} \quad t=0, z \ge 0 : T=T_i \tag{5-2} \end{align*} \] \[ \begin{align*} \left. \eqalign{ {\rm B.C.} \quad t>0,z=0 : & \quad T=T_w \cr z=\infty : & \quad T=T_i \cr } \right\} \tag{5-3} \end{align*} \] ここで、\(T_i\)は固体の初期温度、\(T_w\)は固体の表面温度で、ともに\(t,z\)には依存しない定数とする。

先頭に戻る

例題５の解法その１

　まず解く前に、次の無次元数\(\theta\)を定義し、基礎式、初期条件、境界条件を変換と、 \[ \begin{align*} \theta & = \frac{T-T_i}{T_w-T_i} \tag{5-4} \cr \frac{\partial \theta}{\partial t} & =\lambda \frac{\partial^2 \theta}{\partial z^2} \tag{5-5} \end{align*} \] \[ \begin{align*} {\rm I.C.} \quad t=0, z \ge 0 : \quad \theta=0 \tag{5-6} \end{align*} \] \[ \begin{align*} \left. \eqalign{ {\rm B.C.} \quad t>0,z=0 : & \quad \theta=1 \cr z=\infty : & \quad \theta=0 \cr } \right\} \tag{5-7} \end{align*} \] となる。

　解析解を得るために、次のBoltzmann's transformation（ボルツマン変換）3),4)と呼ばれる変数\(\eta\)を導入して解く方法がある。 \[ \begin{align*} \eta = \frac{z}{2}/{\sqrt {\lambda t}} \tag{5-8} \end{align*} \] 　(5-5)式の両辺の偏微分を\(\eta\)で表すと、 \[ \begin{align*} \frac{\partial \theta}{\partial t} & = -\frac{z}{4\sqrt{\lambda t}} \frac{d \theta}{d \eta} \cr \frac{\partial \theta}{\partial z} & = -\frac{1}{2\sqrt{\lambda t}} \frac{d \theta}{d \eta} \cr \frac{\partial^2 \theta}{\partial z^2} & = \frac{\partial}{\partial z} \Bigl( \frac{\partial \theta}{\partial z} \Bigr) = \frac{1}{4 \lambda t} \frac{d^2 \theta}{d \eta^2} \end{align*} \] となって、(5-5)式は、 \[ \begin{align*} \frac{d^2 \theta}{d \eta^2} + 2\eta \frac{d \theta}{d\eta} = 0 \tag{5-9} \end{align*} \] と、常微分方程式に変換される。一方、初期条件、境界条件は \[ \begin{align*} {\rm I.C.} \quad \eta= \infty : & \quad \theta=0 \tag{5-10} \end{align*} \] \[ \begin{align*} \left. \eqalign{ {\rm B.C.} \quad \eta=0 : & \quad \theta=1 \cr \eta=\infty : & \quad \theta=0 \cr } \right\} \tag{5-11} \end{align*} \] となる。ここで、初期条件と２つ目の境界条件は一致する。(5-9)式を解くと、一般解は \[ \begin{align*} \theta= C_1 \int_0^\eta \exp(-\eta^2)d\eta +C_2 \tag{5-12} \end{align*} \] となる。\(C_1,C_2\)は積分定数。これに初期条件、境界条件を代入する。 \[ \begin{align*} 1 & = C_1 \int_0^0 \exp(-\eta^2)d\eta +C_2 = 0 + C_2 \cr 0 & = C_1 \int_0^\infty \exp(-\eta^2)d\eta +C_2 \cr \end{align*} \] \( \int_0^\infty \exp(-x^2)dx=\sqrt \pi/2\)を考えると、積分定数はつぎのように求まる。 \[ \begin{align*} C_1 &= -\frac{2}{\sqrt \pi} \cr C_2 &= 1 \cr \end{align*} \] したがって、解は \[ \begin{align*} \theta & = 1- \frac{2}{\sqrt \pi} \int_0^\eta \exp(-\eta^2)d\eta \tag{5-13} \cr & = 1 - {\rm erf}\Bigl(\frac{z}{2 \sqrt {\lambda t}}\Bigr) \end{align*} \] となる。無次元化温度を元に戻せば、 \[ \begin{align*} \theta = \frac{T-T_i}{T_w-T_i}=1-{\rm erf}\Bigl(\frac{z}{2\sqrt{\lambda t}}\Bigr) \tag{5-14} \end{align*} \] となる。ここで、\({\rm erf}(x)\)はガウスの誤差関数(error function)と呼ばれ、次式で定義される関数である。また、 \({\rm erfc}(x)\)は余誤差関数(complementary error function)と呼ばれ、次式のように類似の関数で定義される。 \[ \begin{align*} {\rm erf} x & =\frac{2}{\sqrt{\pi}} \int ^x_0 e^{-\xi^2} d \xi \tag{5-15} \cr {\rm erfc} x & =1-{\rm erf}x=\frac{2}{\sqrt{\pi}} \int ^\infty_x e^{-\xi^2} d \xi \tag{5-16} \end{align*} \] 誤差関数には次の性質がある。 \[ \begin{align*} {\rm erf}(0)=0,\quad {\rm erf}( \infty)=1, \quad {\rm erf}(-x)=-{\rm erf}(x) \tag{5-17} \end{align*} \]

　ラプラス変換を紹介していて、別の解法を述べたのは解析解に現われた誤差関数などを紹介したかったからです。 また、円筒座標系のときの解析解にベッセル関数などの数学関数はよく出てきます。
これら関数は、MS-Excelのシート上 からも簡単に呼び出すことができます。ヘルプの「エンジニアリング関数」などで検索できます。

先頭に戻る

例題５の解法その２

　先に示すBoltzmann変換による方法でなく、ラプラス変換を利用した解法を以下に示す。
(5-5)から(5-7)式までは上と同じとする。いま、\( \theta(z,t) \)の時間領域をつぎのラプラス変換 \( \mathcal{L} \theta(z,t) = \Theta(z,s) \)とすると、(5-5)式の左辺は \[ \begin{align*} \frac{\partial \theta}{\partial t}=s \Theta -\theta(z,0)=s\Theta \tag{5-18} \end{align*} \] と変換される。ここで初期条件(5-6)式を用いた。(5-5)式は次の常微分方程式に変換される。 \[ \begin{align*} s\Theta = \lambda \frac{d^2 \Theta}{d z^2} \cr \frac{d^2 \Theta}{d z^2}-\frac{s}{\lambda} \Theta=0 \cr \frac{d^2 \Theta}{d z^2}-q^2 \Theta=0 \tag{5-19} \cr \end{align*} \] 最後の式で、\( q^2=s/\lambda \)とした。境界条件もラプラス変換すると、 \[ \begin{align*} \left. \eqalign{ {\rm B.C.} \quad z=0 : & \quad \Theta=1/s \cr z=\infty : & \quad \Theta=0 \cr } \right\} \tag{5-20} \end{align*} \] (5-19)式の一般解は、 \[ \begin{align*} \Theta = C_1e^{qz}+C_2e^{-qz} \tag{5-21} \end{align*} \] となる。(5-20)式の境界条件を当てはめると \[ \begin{align*} \frac{1}{s} & = C_1+C_2 \cr 0 & = C_1e^{\infty}+C_2e^{-\infty} \end{align*} \] したがって、\(C_1=0,\ C_2=1/s \)が得られる。 \[ \begin{align*} \Theta = \frac{1}{s}e^{-qz} \tag{5-22} \end{align*} \] これをラプラス逆変換することで、求める解析解が得られる。 \[ \begin{align*} \mathcal{L}^{-1}\Theta & = \mathcal{L}^{-1}\frac{e^{-qz}}{s} \cr & ={\rm erfc} \frac{z}{2\sqrt{\lambda t}} \cr & =1- {\rm erf} \frac{z}{2\sqrt{\lambda t}} \tag{5-23} \cr \end{align*} \] 最後の逆変換は、図書4)の変換表を参考にした。
　先のBoltzmann変換に比べれば、はるかにラプラス変換の方がスマートだと思うけど、変換表を準備しないといけない。

先頭に戻る

Literature Cited

　本ページを作成するにあたって、以下の書籍を参考にした。

先頭に戻る