熱流体力学って難しすぎない?

こんにちは，手短に最近あった出来事を紹介していきます．

※これから説明する熱流体力学において初心者であるため，間違った内容が掲載されている可能性があります

・概要

とある授業で

非定常熱伝導方程式

を解かなければならない場面がありました．

ということで概要を話していきます．
授業ではパリレンコーティングした温度センサの応答性について調べました．
また，パリレンの温度特性を調べるために反無限物体の非定常温度場における熱伝導方程式を解きました．

・参考資料

この参考資料では

要点がまとめられている

と感じました．
用語や数式の知識がなく，読むのが少々辛かったのは内緒...

・報告書

今回のまとめを報告書として記載しました．
結果として，

パリレンの熱伝導は温度センサの応答時間に影響を与えない

ことがわかりました．

・プログラム(heat_eq_unsteady.c)

リスト1として図2のグラフのデータを出力するプログラムheat_eq_unsteady.cを記載しました．

//リスト1   heat_eq_unsteady.c
#include<stdio.h>
#include<math.h>

double Simpson(double,double,int); //simpson法によって積分を計算する関数
double IntegralFunction(double); //積分する関数を戻り値として返す関数

int main(){
  //お湯の温度T0=60+273.15[K],パリレンあり温度センサ(室温)T1=...+273.15[K]
  //パリレン内部の温度Ti,時間(変数)t,ガウスの誤差関数の積分上限範囲y=x/(2*sqrt(a*t))
  //定数a=k/(p*c)(ただし、k=熱伝導率,p=密度,c=比熱),j=xとするため、y,aはコメントアウト
  //ガウスの誤差関数erf
  double T0,T1,Ti,t,erf[201];
  //double y,a;
  int kizami = 100,k = 0,j;
  T0=60+273.15;
  T1=25+273.15;
  
  //結果を出力
  for(j=0;j<=200;j++){
    erf[j] = 2 /sqrt(M_PI) * Simpson(0.0,(double)j/100,kizami);
    printf("%.2f,%.3f\n",(double)j/100,erf[j]);
  }

  return 0;
}

//simpson法によって積分を計算する関数
double Simpson(double a,double b,int N){
  double ss1 = 0.0,ss2 = 0.0,x;
  int i;
  for(i=1;i<=N/2-1;i++){
    x = a + (b-a) / (double)N * (double)i*2.0;
    ss1 += IntegralFunction(x);
  }
  for(i=1;i<=N/2;i++){
    x = a + (b-a) / (double)N * (double)i*2.0;
    ss2 += IntegralFunction(x);
  }
  return (b-a)/(3.0*(double)N) * (IntegralFunction(a) + IntegralFunction(b) + 2.0 * ss1 + 4.0 * ss2);
}

//積分する関数を戻り値として返す関数
double IntegralFunction(double a){
  return exp(-pow(a,2));  
}

・プログラム(TheoreticalValue_temp.c)

リスト2として図3の温度応答を計算するプログラムTheoreticalValue_temp.cを記載しました．
プログラム名のイニシャルが大文字になっている...

/*
リスト2   TheoreticalValue_temp.c
図2のグラフのデータを出力するプログラムheat_eq_unsteady.cで出力されたデータheat_eq_unsteady.txtを
もとにパリレンCの温度特性を調べる．
*/

#include<stdio.h>
#include<math.h>
#include<stdlib.h>

int main(){
  int i=0,j,k,dT;
  //お湯の温度T0=80.0度,初期温度Ti=25.0度,定数a=,パリレン幅d=25μm
  //定数a=K/(pc),熱伝導率K=0.084W/(m*K),比熱c=712J/(kg11*K),密度p=1289kg/m^3
  double T0=80.0,Ti=25.0,K=0.084,c=712,p=1289,d=25e-6;
  double num,min,T,x,y,time,a,n1,n2;
  double data1[201],data2[201],sub[201];
  FILE *fp;

  fp = fopen("heat_eq_unsteady.txt","r");

  if(fp == NULL){
    printf("ファイルを開けませんでした\n");
    return -1;
  }

    
  //温度の入力
  printf("任意の温度T(x,u)を入力して下さい\n");
  scanf("%lf",&T);
  

  //データを読み取る
  while((fscanf(fp,"%lf,%lf",&n1,&n2)) != EOF){
    data1[i] = n1;
    data2[i] = n2;
    i++;
  }
  
   //入力した温度に対するグラフの座標を求める                                                                                                                                                              
    y = (T-T0)/(Ti-T0);
    for(j=0;j<i;j++){
      if(y > data2[j]){
        sub[j] = y - data2[j];
      }
      else{
        sub[j] = data2[j] - y;
      }
    }
    min = sub[0];
    for(j=1;j<i;j++){
      if(sub[j]<min){
        min = sub[j];
        k = j;//data2[k]がyともっとも近い値となる                                                                                                                                                           
      }
    }
    //y=data2[k]に対応したxを求める                                                                                                                                                                         
    x = data1[k];

    a = K/(p*c);
	
	//時間を計算し，結果として出力
    time = pow(d,2)/(4*a*pow(x,2));
    printf("%.3f\n",time);
   
  fclose(fp);
  return 0;
}