右上の添字は時間を、右下の添字は空間を表すのが慣例なので、そのようにし、時間項をオイラー陽解法で、2階微分を2次中心差分で離散化します。\label{eq:scale-factor-2}\end{align}$$境界条件と初期条件を考えます。\label{eq:scale-factor-4}\end{align}$$となります。ここまではいいでしょう。このように時間を前進差分で離散化し、空間を中心差分で離散化した場合には、次のn+1時点の値を下図のように、n時点での両隣の値と自身の値から求めることになります。空間の端っこの値をえいやっと決めなければ、解くことはできません。いった意味で空間の端っこの値を決める条件を境界条件と言っています。最初の値を空間全体で決めてやらなければ、時間的に次の値を求めることができないので、解くことはできません。space法)と呼ぶらしいです。一応行列で表現するとどうなるか、書いてみます。\right]\end{align}$$できました!あとは、この連立方程式を解けば終わりです。今回の連立方程式はまじでそのまま解けます。素直にやっていってあとは繰り返すだけです。トーマス法でやった時との違いがわかるでしょうか?よくみてください。トーマス法の時に使った行列がこれです。\right]\end{align}$$求めたいものを、\(\mathbf{Y}\)、すでにわかっているものを\(\mathbf{X}\)としましょう。\mathbf{X}\)となっており、逆行列を解くような作業をしていました。トーマス法は特殊条件の時のみ厳密に逆行列が解ける方法です。一般的には、いろんな解法で逆行列を解かねばなりません。\mathbf{Y}\)という形になってます。プログラムではプログラムを実装しましょう。言語はとりあえずCでやります。他の言語も暇だったらやります。以下のプログラムは、冒頭でも紹介した「流体解析の基礎」からデバッグ・改変して使っています。元々のプログラムだと、アウトプットもなくあまり解いてる感もないので、アニメーション作ります。\end{eqnarray}という形をとりあえず使います。以下のファイルをHEAT_E.cppとでもして保存しましょう。ここでパラメータを設定します。そのそも時間も長さも単位を定義していないので、単位はあまり考えない。sys.stat.hがないと機能しない。その場合は手動でanimeのフォルダを作ってコメントアウトする。ここでパラメータを設定します。そのそも時間も長さも単位を定義していないので、単位はあまり考えない。sys.stat.hがないと機能しない。その場合は手動でanimeのフォルダを作ってコメントアウトする。"\t%d\t%f\n",i,uw[i]);}fclose(file);}解説です。解説ファイルがついているので、それを読みます。解くプログラムである。プログラムは非常に単純であり、\eqref{eq:scale-factor-4}式をそのまま記述すればよい。uはxとtの関数であるため二次元の配列をとってもよいが、時間発展型の問題であり途中の時間ステップでのuの値を残しておく必要がない場合には一次元配列で十分である。(uuとする)を記憶する配列を二種類用意して\eqref{eq:scale-factor-4}でuからuuを各格子点で計算したあと、uuをuでおきかえる。プログラムのはじめの部分において計算に必要なパラメータを読んだり、計算したりする。熱伝導率は1にしているが、1でない場合にはrの右辺に熱伝導率をかければよい。i5は表示のために用いる。刻みに結果を表示する。CONDITIONのループで初期条件を与えている。1)を初期条件に用いているが、別の関数を用いる場合にはこの部分を書き換える。LOOPのところのメインループで熱伝導方程式を解いている。メインループの行のすぐ下の2行が境界条件である。境界において関数値は固定されているから一度指定すれば十分であり、ループの中に入れる必要はないが、条件によっては時間的に変化する場合もあるのでこの部分で指定している。おきに、出力のためのサブルーチンを呼んだあと、次のループで\eqref{eq:scale-factor-4}の右辺を計算し配列uuに代入し、その次のループでuuをuに置き換えてひとつの時間ステップを終了する。今回はついでに、僕が作ったアニメーション用の書き出し関数output_animeでアニメーションを書きます。格子点のX座標らしいです。コード置いときますこの微分方程式を解いて熱の拡散を非定常に、つまり0.005とか進むたびに全体の温度を求めているということになります。計算を回してみましょう。HEAT_E.cpp./a.outで実行できるかと思います。デフォルトだとa.outになるので、名前をつけたい場合は-oオプションでできます。heat_e.out./heat_e.out実行すると、animeというフォルダが作成されて、そこにファイルが書き出されると思います。これをgnuplotでアニメーションにしていきます。めんどいので、gnuplotのプログラムを作るプログラムで作ります。以下のプログラムをheat_animetion_make.cppとでもして保存してください。formatを使うためのtemplate。メッシュ数に合わせる。formatを使うためのtemplate。メッシュ数に合わせる。1"<<endl;}}file.close();}こいつで、anime.gpといった奴を作り、それを使ってgnuplotで簡単にアニメーションを書きます。anime.gpだんだん、熱が逃げていく形になっているのがわかるでしょうか?今回gnuplotのアニメーションの描画雑めです。ループ使って直で書けるはずなんですが、意外とむずいんでこんなんで。あとは、初期条件とか境界条件いじってみます。x;}}このあたりをいじって境界条件を変えたりできます。初期条件はまあいいので、境界条件を変えます。\label{eq:scale-factor-5}\end{align}$$とします。0\)のように、勾配値で与える境界条件をノイマン境界条件といいます。それぞれの境界条件の意味を考えてみます。\(u(0)=0\)は固定条件です。今回は熱ですから、熱が0度の物体or空気に接している、といった状況が想定できます。言うと、内部からの熱の影響は受けるはずなので、それを受けないくらいの規模感のものに接していると考えられます。0\)の勾配条件はどうでしょうか?端っこで勾配がゼロと言っているので、隣の値とだいたい同じくらいになりそうではないでしょうか?隣の値に引きづられて、同じ値になっていくような挙動になります。1\)と言ったような条件だとどうなるでしょうか?右端での傾きがプラスとなり、つまりは左の一個隣の点よりも高い温度となり、左の一個隣の点に対して、常に熱を与えているような条件になります。これは、flux(熱流束)の条件であり、ヒーター等で熱を与え続けているイメージです。左だと-1が与える条件になり、ゴリゴリに手で書く時は左右で違うので面倒です。これらに対し、プログラムでどのように表現すればいいかを考えます。右端での、境界を考え、\(u_J\)が端っこの点の値、\(u_P\)は領域外の仮想的な点の値です。今回は中心差分で考えていたので、境界条件も中心差分で考えます。\label{eq:scale-factor-6}\end{align}$$となります。\(u_J\)での差分方程式も考えます。\label{eq:scale-factor-7}\end{align}$$ここから、\(u_P\)を消去します。\label{eq:scale-factor-8}\end{align}$$これが、今回適用する境界条件になります。2*r)*u[mx];こうします。1\)も考えてみます。\label{eq:scale-factor-10}\end{align}$$となります。最後の項の分が常に与え続けられることになります。1\)の単位が[K/m]となっています。1[K/m]を与えていますが、これに流体の熱伝達率[W/m・K]をかけてやれば、熱流速の単位[W/m2]になります。今回は1*(熱伝達率)[W/m2]の熱量を与えていたことになります。1次元だから、ちょっと違うかもしれませんが、まあ目をつむってください。2*r*dx;に書き換えれば終わりです。微妙な疑問として、uu[mx]に対して書き換えるのもありなんじゃないか、つまり境界条件を先に変える?後に変える?というのもあるのですが、境界条件を先に変えるのが一般的なようです。そこまで違いは出ないと思いますが。流体を学ぶ前に、流体とは何かから行きましょう。プログラムは、以下の書籍の内容やプログラムを参考にして作ってます。コメントをキャンセルメールアドレスが公開されることはありません。