1次元熱伝導を陽解法で解く 一次元の熱伝導を解きます

非圧縮性流体

まずは陽解法で

プログラムは、以下の書籍の内容やプログラムを参考にして作ってます。プログラムは、一部改修・デバッグして使用しています。超いい本なので、ぜひ読んでください。

プログラムは朝倉書店のホームベーじからダウンロードできます。書籍名で検索してみてください。

コメント書いたり修正したプログラムも、githubにアップしておきます。

 

今回は、時間項が入った方程式で要は非定常解析の一番簡単なもの、1次元熱伝導方程式を解きます。

$$
\begin{align}
&\frac{ \partial u }{ \partial t}  = k\frac{\partial^2 u}{\partial x^2} \qquad (k > 0, t > 0,0<x<1) \\\\
& u(0,t) = u(1,t) = 0 \\\\
& u(x,0) = f(x)
\tag{1} \label{eq:scale-factor-1}
\end{align}
$$

を考えます。

前回と同じように、以下のようにメッシュを切ります。今回は、時間に関しても離散化します。

J-1個の格子、J個の格子点に離散化します。格子点の番号を左から1,…,Jとし、n時点での、\(u\)の値を、\(u_{1}^{n},u_{2}^{n},..,\)、と名前をつけていきます。代表点を、\(u_{j}^{n}\)とします。

前回の格子点は0から数えましたが、今回は1から数えています。ややこしくなってすみません。ただ、どちらの流派もあるので、慣れておくといいと思います。時間は制限はありません。どこまでもやりたいとこまでやります。

なので、番号はx=0の点で格子番号は1, x=1の点で格子番号はJ, t=0の点で時間の方の格子番号は1になってます。ややこしいです。

格子幅を等間隔で\(\Delta x, \Delta t \)とします。右上の添字は時間を、右下の添字は空間を表すのが慣例なので、そのようにし、時間項をオイラー陽解法で、2階微分を2次中心差分で離散化します。

$$
\begin{align}
&\frac{u_{j}^{n+1} – u_{j}^{n}}{\Delta t} = k \frac{u_{j-1}^{n} – 2u_{j}^{n} + u_{j+1}^n}{(\Delta x)^2} \\\\
& u_{j}^{n+1} = r u_{j-1}^{n} + (1-2r) u_{j}^{n} +r u_{j+1}^n \qquad (r = k \Delta t/ (\Delta x)^2) \\\\
&(j=2,,…,J-1, n=1,2,…)
\tag{2} \label{eq:scale-factor-2}
\end{align}
$$

境界条件と初期条件を考えます。

$$
\begin{align}
& u(0,t) = u(1,t) = 0 \\\\
& u(x,0) = f(x) \\\\\
& より、\\\\
& u_{1}^{n} = u_{J}^n = 0 \qquad (n=1,2,…)\\\\
& u_{j}^{1} = f(x_j) \qquad (j=1,2,…,J)
\tag{4} \label{eq:scale-factor-4}
\end{align}
$$

となります。ここまではいいでしょう。

このように時間を前進差分で離散化し、空間を中心差分で離散化した場合には、次のn+1時点の値を下図のように、n時点での両隣の値と自身の値から求めることになります。なので空間の端っこの値をえいやっと決めなければ、解くことはできません。こういった意味で空間の端っこの値を決める条件を境界条件と言っています。

また、一番最初の値を空間全体で決めてやらなければ、時間的に次の値を求めることができないので、解くことはできません。時間的な最初の条件を初期条件と言っています。

時間前進差分、空間中心差分の方法を、FTCS(forword time center space法)と呼ぶらしいです。

一応行列で表現するとどうなるか、書いてみます。

$$\begin{align}
& \left[ \begin{array}{cccccc} 1-2r & r &  &  &  & 0 \\ r & 1-2r & r &  &  & \\  & \ddots & \ddots & \ddots &  &  \\  &  & r & 1-2r & r & \\ 0 &  &  &  & r & 1-2r \end{array} \right] \hspace{ 3pt } \left[ \begin{array}{c} u_{2}^{n} \\ u_{3}^{n} \\ \vdots \\ u_{J-2}^{n} \\ u_{J-1}^{n} \end{array} \right] \hspace{ 3pt }=\hspace{ 3pt } \left[ \begin{array}{c} u_{2}^{n+1} \\ u_{3}^{n+1} \\ \vdots \\ u_{J-2}^{n+1} \\ u_{J-1}^{n+1} \end{array} \right]
\end{align}$$

できました!あとは、この連立方程式を解けば終わりです。


解法

今回の連立方程式はまじでそのまま解けます。素直にやっていってあとは繰り返すだけです。

トーマス法でやった時との違いがわかるでしょうか?よくみてください。トーマス法の時に使った行列がこれです。

$$\begin{align}
& \left[ \begin{array}{cccccc} \Delta x^2 -2 & 1 &  &  &  & 0 \\ 1 & \Delta x^2 -2 & 1 &  &  & \\  & \ddots & \ddots & \ddots &  &  \\  &  & 1 & \Delta x^2-2 & 1 & \\ 0 &  &  &  & 1 & \Delta x^2 -2 \end{array} \right] \hspace{ 3pt } \left[ \begin{array}{c} u_{1} \\ u_{2} \\ \vdots \\ u_{N-2} \\ u_{N-1} \end{array} \right] \hspace{ 3pt }=\hspace{ 3pt } \left[ \begin{array}{c} -\Delta x^3 \\ – 2 \Delta x^3 \\ \vdots \\ -(N-2)\Delta x^3 \\ -(N-1)\Delta x^3 \end{array} \right]
\end{align}$$

求めたいものを、\(\mathbf{Y}\)、すでにわかっているものを\(\mathbf{X}\)としましょう。

トーマス法の法では、\(A\mathbf{Y} = \mathbf{X}\)となっており、逆行列を解くような作業をしていました。トーマス法は特殊条件の時のみ厳密に逆行列が解ける方法です。一般的には、いろんな解法で逆行列を解かねばなりません。

しかし、今回は、\(u_{j}^{n+1}\)が求めたいものですので、既に、\(A\mathbf{X} = \mathbf{Y}\)

という形になってます。なら、やるだけですね。


プログラム

ではプログラムを実装しましょう。言語はとりあえずCでやります。他の言語も暇だったらやります。

以下のプログラムは、冒頭でも紹介した「流体解析の基礎」からデバッグ・改変して使っています。

元々のプログラムだと、アウトプットもなくあまり解いてる感もないので、アニメーション作ります。

そして、初期条件に、

\begin{eqnarray} f(x)  = \begin{cases} x & (0 \leq x \leq 0.5 ) \\ 1-x & (0.5 \leq x \leq 1 ) \end{cases} \end{eqnarray}

という形をとりあえず使います。

以下のファイルをHEAT_E.cppとでもして保存しましょう。

 

解説です。まず、解説ファイルがついているので、それを読みます。

<一次元熱伝導方程式(FTCS法)>

 このプログラムはFTCS法で一次元熱伝導方程式を解くプログラムである。プログラムは非常に単純であり、\eqref{eq:scale-factor-4}式をそのまま記述すればよい。uはxとtの関数であるため二次元の配列をとってもよいが、時間発展型の問題であり途中の時間ステップでのuの値を残しておく必要がない場合には一次元配列で十分である。なお一次元配列を用いる場合には現時間ステップでの\(u^n\) (uとする)と次の時間ステップでの \(u^{n+1}\) (uuとする)を記憶する配列を二種類用意して\eqref{eq:scale-factor-4}でuからuuを各格子点で計算したあと、uuをuでおきかえる。
 プログラムのはじめの部分において計算に必要なパラメータを読んだり、計算したりする。mxはx方向の格子数、kmは時間ステップ数、dtは時間刻み\(Delta t\)である。なお熱伝導率は1にしているが、1でない場合にはrの右辺に熱伝導率をかければよい。i5は表示のために用いる。(t=0.05 刻みに結果を表示する。)最初INITIAL CONDITIONのループで初期条件を与えている。ここでは
  f(x) = x    (0 ≦x ≦ 0.5)
       = 1- x  (0.5 <x ≦ 1)
を初期条件に用いているが、別の関数を用いる場合にはこの部分を書き換える。MAIN LOOPのところのメインループで熱伝導方程式を解いている。メインループの行のすぐ下の2行が境界条件である。境界において関数値は固定されているから一度指定すれば十分であり、ループの中に入れる必要はないが、条件によっては時間的に変化する場合もあるのでこの部分で指定している。無次元時間でkop おきに、出力のためのサブルーチンを呼んだあと、次のループで\eqref{eq:scale-factor-4}の右辺を計算し配列uuに代入し、その次のループでuuをuに置き換えてひとつの時間ステップを終了する。

今回はついでに、僕が作ったアニメーション用の書き出し関数output_animeでアニメーションを書きます。

メインプログラムの仕様

  配列 u —- 現時点での解を記憶する配列
  配列 uu —- 前の時間ステップでの解を記憶する配列
  変数 nx —- プログラムの配列の大きさを指定するパラメータ
  変数 kx —- 格子数
  変数 mx —- 格子点数
  変数 km —- 計算打ち切りのタイムステップ数
  変数 dt —- 時間間隔
  変数 dx —- 格子間隔
  変数 r —- (時間間隔)/(格子間隔)**2
変数 ih —- 中央の格子点番号
  変数 i5 —- 無次元時間0.05に対応するタイムステップ数
  変数 x —- 格子点のX座標

らしいです。

コード置いときます

この微分方程式を解いて熱の拡散を非定常に、つまり0.005とか進むたびに全体の温度を求めているということになります。

まずは、計算を回してみましょう。

で実行できるかと思います。デフォルトだとa.outになるので、名前をつけたい場合は-oオプションでできます。

実行すると、animeというフォルダが作成されて、そこにファイルが書き出されると思います。

これをgnuplotでアニメーションにしていきます。

ただ、色々とめんどいので、gnuplotのプログラムを作るプログラムで作ります。

以下のプログラムをheat_animetion_make.cppとでもして保存してください。

 

こいつで、anime.gpといった奴を作り、それを使ってgnuplotで簡単にアニメーションを書きます。

だんだん、熱が逃げていく形になっているのがわかるでしょうか?

今回gnuplotのアニメーションの描画雑めです。一応gnuplotはfor ループ使って直で書けるはずなんですが、意外とむずいんでこんなんで。

あとは、初期条件とか境界条件いじってみます。

このあたりをいじって最初のxとかihとか1 -xの部分をいじったり、

このあたりをいじって境界条件を変えたりできます。

 


境界条件の変更

初期条件はまあいいので、境界条件を変えます。

$$
\begin{align}
& u(0) = 0\\\\
& \left.\frac{\partial u}{\partial x}\right\vert_{x=1} = 0 \\\\
\tag{5} \label{eq:scale-factor-5}
\end{align}
$$

とします。境界条件には大別して2種類あり、\(u(0)=0\)のように、固定値で与える境界条件をディリクレ境界条件といい、\(\left.\frac{\partial u}{\partial x}\right\vert_{x=1} = 0\)のように、勾配値で与える境界条件をノイマン境界条件といいます。

それぞれの境界条件の意味を考えてみます。

まず、\(u(0)=0\)は固定条件です。今回は熱ですから、熱が0度の物体or空気に接している、といった状況が想定できます。もっと言うと、内部からの熱の影響は受けるはずなので、それを受けないくらいの規模感のものに接していると考えられます。

\(\left.\frac{\partial u}{\partial x}\right\vert_{x=1} = 0\)の勾配条件はどうでしょうか?端っこで勾配がゼロと言っているので、隣の値とだいたい同じくらいになりそうではないでしょうか?隣の値に引きづられて、同じ値になっていくような挙動になります。

\(\left.\frac{\partial u}{\partial x}\right\vert_{x=1} = 1\)と言ったような条件だとどうなるでしょうか?右端での傾きがプラスとなり、つまりは左の一個隣の点よりも高い温度となり、左の一個隣の点に対して、常に熱を与えているような条件になります。これは、flux(熱流束)の条件であり、ヒーター等で熱を与え続けているイメージです。左だと-1が与える条件になり、ゴリゴリに手で書く時は左右で違うので面倒です。

これらに対し、プログラムでどのように表現すればいいかを考えます。

まずは、\(\left.\frac{\partial u}{\partial x}\right\vert_{x=1} = 0\)の条件から。

右端での、境界を考え、\(u_J\)が端っこの点の値、\(u_P\)は領域外の仮想的な点の値です。

今回は中心差分で考えていたので、境界条件も中心差分で考えます。すると、

$$
\begin{align}
& \frac{u_P – u_{J-1} }{2 \Delta x} = 0 \\\\
\tag{6} \label{eq:scale-factor-6}
\end{align}
$$

となります。また、\(u_J\)での差分方程式も考えます。

$$
\begin{align}
& u_{J}^{n+1} = ru_{J-1}^{n} + (1-2r)u_{J}^{n} + ru_{P}^{n} \\\\
\tag{7} \label{eq:scale-factor-7}
\end{align}
$$

ここから、\(u_P\)を消去します。\eqref{eq:scale-factor-7}式から、\(u_P = u_{J-1}\)となるので、

$$
\begin{align}
& u_{J}^{n+1} = 2ru_{J-1}^{n} + (1-2r)u_{J}^{n}
\tag{8} \label{eq:scale-factor-8}
\end{align}
$$

これが、今回適用する境界条件になります。

プログラム上では、

この部分を、

こうします。

 

もう一個、\(\left.\frac{\partial u}{\partial x}\right\vert_{x=1} = 1\)も考えてみます。

同じように、

$$
\begin{align}
& \frac{u_P-u_{J-1}}{2 \Delta x} = 1 \\\\
\tag{9} \label{eq:scale-factor-9}
\end{align}
$$

より、\(u_P = 2\Delta x + u_{J-1}\)となるので、

$$
\begin{align}
& u_{J}^{n+1} = 2ru_{J-1}^{n} + (1-2r)u_{J}^{n} + 2 r \Delta x  \\\\
& u_{J}^{n+1} = 2ru_{J-1}^{n} + (1-2r)u_{J}^{n} + 2 k \frac{\Delta t}{\Delta x}  \\\\
\tag{10} \label{eq:scale-factor-10}
\end{align}
$$

となります。最後の項の分が常に与え続けられることになります。単位を考えてみると、最初に与えた、\(\left.\frac{\partial u}{\partial x}\right\vert_{x=1} = 1\)の単位が[K/m]となっています。1[K/m]を与えていますが、これに流体の熱伝達率[W/m・K]をかけてやれば、熱流速の単位[W/㎡]になります。

つまり、今回は1*(熱伝達率)[W/㎡]の熱量を与えていたことになります。1次元だから、ちょっと違うかもしれませんが、まあ目をつむってください。

に書き換えれば終わりです。

微妙な疑問として、uu[mx]に対して書き換えるのもありなんじゃないか、つまり境界条件を先に変える?後に変える?というのもあるのですが、境界条件を先に変えるのが一般的なようです。そこまで違いは出ないと思いますが。



[
前の記事へ]  [非圧縮性流体の目次へ]  [次の記事へ]