#はじめに
初投稿です。ふと数値計算をしたくなったので熱伝導シミュレーションのプログラミングを作ってみました。数値計算のスクリプトにはC言語、グラフの作成にはgnuplotを用いました。
####作るもの
今回目標としたのは2次元平面で任意の場所に任意の個数熱源を置き、その熱伝導の時間発展をプロットすることを目的としました。利用した差分法はFTCSです。
初投稿で至らぬ部分や足りない部分等ありますが随時加筆していきます。
#方針
少し長くなるので最初に目次を立てます。
- 準備
- 数理モデルを考える
- Cによるコードの記述
- gnuplotによるグラフのプロット
#準備
今回はC言語で数値計算を行い、その結果を.csvファイルに保存して、gnuplotでプロットすることを目標としています。よって必要なのはプログラムを書くためのエディター(今回はVScodeを使いました。)とgnuplotのインストールです。
gnuplotのインストールはこちらの方の記事が参考になります。
※すいませんまだ参考になる記事見つけられていません
#数理モデル
###FTCSとは
FTCSとはFoward Time Centered Spaceの略であり、次のような方法を用います。
まず、以下の偏微分方程式を考えます。
$$
\frac{\partial u}{\partial t} = F(u, x, t , \frac{\partial^2 u}{\partial x^2})
$$
ここで、$u(x,t) = u(i \Delta x, n\Delta t)= u^n_i$と置くと、左辺を微分の定義に従って展開することで
$$
\frac{u^{n+1}_i - u^n_i}{\Delta t} = F^n_i(u, x, t,\frac{\partial^2 u}{\partial x^2} )
$$
を得られます。この操作は右辺で前進差分をとったということになります。
さて、ここで次の式で与えられる一次元熱拡散方程式を考えます。
$$
\frac{\partial u}{\partial t} = \alpha \frac{\partial^2 u}{\partial x^2}
$$
右辺について考えると、$u(x \pm \Delta x,t)$を$\Delta x$が2次の項まで考えることで
$$
\frac{\partial^2 u}{\partial x^2} \simeq
\frac{u^n_{i+1} - 2u^n_i + u^n_{i-1}}{\Delta x ^2}
$$
を得られます。この形はtの偏微分を分解したときに得た前進差分とは違い、$i+1,i,i-1$の項を用いているので中央差分と呼ばれます。
このようにして時間に対しては前進差分、空間に対しては中央差分をとった式は
$$
\frac{u^{n+1}_i - u^n_i}{\Delta t} =
\alpha \frac{u^n_{i+1} - 2u^n_i + u^n_{i-1}}{\Delta x ^2}
$$
となり、このような差分のとり方はFTCSスキームと呼ばれます。
###2次元熱拡散方程式
ここでは $u(x,y,t) = u(i\Delta x, j\Delta y, n\Delta t) = u^n_{i,j}$ とし、一次元の熱拡散方程式を二次元に拡張します。拡散係数を$D$と置くと結果的に式は
$$
u^{n+1}_{i,j} = D(\frac{u^n_{i+1,j} - 2u^n_{i,j} + u^n_{i-1,j}}{\Delta x ^2} + \frac{u^n_{i,j+1} - 2u^n_{i,j} + u^n_{i,j-1}}{\Delta y ^2})
$$
となります。
###計算手法
FTCSによる離散化が終わったので、次は計算手法について考えます。
一番最初のgifのようにあらわしたい場合、各メッシュによる計算結果を座標とともに出力、すなわち座標（２次元）＋温度（1次元）の三次元の情報を各ステップ（時間）で書きださなければいけません。よってメッシュの合計がA個、プロットする回数がN回であった場合、中身が3Aの数値を持つ.csvファイルがN個必要であることになります。
このことから計算を行う前にまずメッシュの数、実際に計算を行う回数と結果をプロットする回数、さらに初期条件と境界条件が必要であることがわかります。
ここで、この計算を行う場合配列を用いた計算が便利そうだと考えたので、JK行列の各成分を温度として計算を行います。（ここで行列を用いましたが行列の演算はほぼ使いません）
#コードの記述
###数値計算スクリプト
まず、おまじないとしてincludeファイルを呼び出します。そしてのちに定義する関数で使うNRの変数を定義しておきます。
#if defined(__STDC__) || defined(ANSI) || defined(NRANSI) /* ANSI */
#include <stdio.h>
#include <stdlib.h>
#define NR_END 1
#define FREE_ARG char*
次にメッシュの情報を格納するための行列を定義します。C言語の純粋な2次元配列では配列の呼び出し、すなわち計算する過程が非効率的なのでoffset化と呼ばれる手法で簡略化します。また、本来配列の順番は0から始まりますが、行列として使いたいので1から始まるように設定を行います。
void nrerror(char error_text[])
/* Numerical Recipes standard error handler */
{
fprintf(stderr,"Numerical Recipes run-time error...\n");
fprintf(stderr,"%s\n",error_text);
fprintf(stderr,"...now exiting to system...\n");
exit(1);
}
/* double型の2次元配列をオフセット化&ポインタの配列へのポインタ変換 */
double **dmatrix(long nrl, long nrh, long ncl, long nch)
/* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */
{
long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
double **m;
/* allocate pointers to rows */
m=(double **) malloc((size_t)((nrow+NR_END)*sizeof(double*)));
if (!m) nrerror("allocation failure 1 in matrix()");
m += NR_END;
m -= nrl;
/* allocate rows and set pointers to them */
m[nrl]=(double *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(double)));
if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
m[nrl] += NR_END;
m[nrl] -= ncl;
for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
/* return pointer to array of pointers to rows */
return m;
}
#endif /* ANSI */
次に偏微分やその初期条件等で用いる定数を定義します。上から微小温度、微小メッシュ幅(x)、微小メッシュ幅(y)、縦メッシュの数、横メッシュの数、初期の全体の温度、熱源の温度、計算をループさせる回数の上限、csvファイルに書き込むスパン、拡散係数 となっています。
ここでの「csvファイルに書き込むスパン」とは、csvに書き込むタイミングのことを指します。例えば、ここでは２０と設定されていますが、20回に一回csvファイルに書き込み結果を出力するということです。全計算ループ回数は2000で与えられているので、出力されるcsvファイルの数は100個となります。
#define DT 1e-4
#define DX 1e-1
#define DY 1e-1
#define latT 100
#define latY 100
#define TEMP0 25
#define TEMP1 300
#define NEND 2000
#define NSAVE 20
#define D 20
さて、次は各変数の定義を行います。最終的に求めたいのは各メッシュの温度の時間発展なので、行列によってそれを定義します。
偏微分方程式には境界条件が必要ですが、今回は境界を自由端として考えます。自由端とは、境界のメッシュが1メッシュ分内側の温度と同じ値を持つということです。
int main(void){
/* i:行 j:列 n:ループ回数 */
int i, j, n;
/* 温度の宣言 */
double **temp, **ntemp;
temp=dmatrix(0, latT+1, 0, latY+1);/* 境界条件を考えるために余分に一回り大きい行列を考える */
ntemp=dmatrix(0, latT+1, 0, latY+1);
/* 残りの変数の定義 */
char fname[256];
int a=20, b=30;/* 熱源の座標 */
int c=80, d=60;
double kappa = D;/* 拡散係数 */
double dx=DX,dy=DY,dt=DT; /* 各微小量 */
変数の定義を終えたら、次は初期条件の設定を行います。ここで境界も同様に初期温度を代入していきます。
この操作の後に、csvファイルにこの時の全メッシュの情報を書き込みます。
/* 各格子点の温度の初期化 */
for(i = 0; i <= latT+1; i++){
for(j =0; j <= latY+1; j++){
temp[i][j] = TEMP0;
if(i==a && j==b) temp[i][j] = TEMP1;
if(i==c && j==d) temp[i][j] = TEMP1;
//printf("%d, %d, %g\n",i, j, temp[i][j]);
}
}
/* 初期値の書き込み */
FILE *fp;
if((fp = fopen("kekka000.csv","wt")) == NULL) exit(1);/* if の中でfpを定義してかつfopenしている？*/
for(i = 1; i <= latT; i++){
for(j =1; j <= latY; j++){
fprintf(fp, "%g,%g,%g\n", dx*j, dy*i, temp[i][j]); ///dxを１としたときにおかしい値になったので1.0と置いた
}
fprintf(fp, "\n"); /* 改行を入れないとgnuがエラーを吐く */
}
fclose(fp);
さて、次はいよいよ熱源を設置して計算を行います。今回設置する熱源は300[℃]のものであり、大体オーブンと同じくらいの温度です。
for(n=1;n<=NEND;n++){
/* 熱源を設置 */
temp[a][b] = TEMP1;
temp[c][d] = TEMP1;
/* 境界条件(自由端)の更新 各頂点は計算過程に影響しないので無視 */
for(i=1;i<=latT;i++){
temp[i][0] = temp[i][1];
temp[i][latY+1] = temp[i][latY];
}
for(j=1;j<=latY;j++){
temp[0][j] = temp[1][j];
temp[latT+1][j] = temp[latT][j];
}
/* 差分法のFTCSによる計算 ここで、(a,b)の温度は一定であることに注意する。*/
for(i=1;i<=latT;i++){
for(j=1;j<=latY;j++){
if(i==a && j==b) ntemp[i][j] = TEMP1;
else if(i==c && j==d) ntemp[i][j] = TEMP1;
else{ntemp[i][j] = temp[i][j] +/* jがx, iがyに対応していることに注意する */
kappa*(
(temp[i][j+1] -2*temp[i][j] + temp[i][j-1])/(dx * dx) +
(temp[i+1][j] -2*temp[i][j] + temp[i-1][j])/(dy * dy)
) * dt;
}
}
}
/* 温度の更新 */
/* 各格子点の温度の初期化 */
for(i = 0; i <= latT+1; i++){
for(j =0; j <= latY+1; j++){
temp[i,j] = ntemp[i,j];
if(i==a && j==b && (n%NSAVE==0 || n==0)) printf("T(%g,%g) = %g\n",j*dx,i*dy, temp[i][j]);
}
}
/* csvへの書き込み */
if(n%NSAVE==0){
sprintf(fname,"kekka%03d.csv",n / NSAVE); /* fname[] にファイル名を代入 */
if((fp = fopen(fname,"wt")) == NULL) exit(1);
for(i = 1; i <= latT; i++){
for(j =1; j <= latY; j++){
fprintf(fp, "%g,%g,%g\n", dx*j, dy*i, temp[i][j]);
}
fprintf(fp, "\n");
}
fclose(fp);
}
}
return 0;
}
###gnuplotによるプロット
#参考記事
https://en.wikipedia.org/wiki/FTCS_scheme
#引用文献
ニューメリカルレシピ・イン・シー 日本語版―C言語による数値計算のレシピ