ティラー

2014年10月2日
本日の講義及び演習
• 偏微分方程式の偏微分項をコンピュータで扱える
ようにする
– 離散化（差分化）
• テイラー展開の利用
数値シミュレーション
– 1階微分項に対する差分式
– 2階微分項に対する差分式
2014年度 第2回
• 1次元熱伝導方程式に適用して差分式を導出
• Excelを利用した温度変化シミュレーション
永野 （熱流体システム研究室）
[email protected]
重要！
2
熱の伝わり方（伝熱モード）
熱伝導方程式の導出
重要！
熱移動量の表現
1次元（x方向のみ）を考える
熱伝導・・・物質（固体・流体）内の熱移動
dQ [ J ]
：熱移動量
dQ [J/s]  [ W ] ：単位時間当たり
面積dA
T1
q
熱伝達（対流）・・・流体と物体間の熱移動
熱の移動方向は
x方向を正とする
熱移動量dQ
流体
温度差と温度勾配
T2
dT  T2  T1 [ K ] ：温度差
dT
[ K/m ] ：温度勾配
dx
＝単位距離当たりの温度差
輻射（放射）・・・電磁波として移動する熱移動
電磁波
距離
dQ
[ W m 2 ] ：単位時間 単位面積当たり
dA
熱流束（heat flux）
dx
x
3
フーリエの法則
復習
4
復習
温度勾配 と 熱流束 の正負
x
フーリエの法則
熱流束（単位時間 単位面積当たりの熱
移動量）は温度勾配に比例する
面積dA
T1
q  
熱移動量dQ
dx
x
－
q
T2
－
T1
①
q
T2
比例係数λを 熱 伝 導 率
（thermal conductivity）と呼ぶ
T2
距離
dT
dx
T1
+
q
熱がxの正方向に流れるのはT1 > T2の時，
すなわちT2－T1 < 0となる時である．
よって負記号が付く．
5
q
+
熱流束
dT
0
dx
q0
dT
0
dx
q0
dT
0
dx
q0
dT
0
dx
q0
T2
T1
T1
温度勾配
T2
6
1
熱伝導率の単位
復習
• 物質を挟んだ点１と点２の間にある温度差があった
とする．
熱伝導率が大きいと移動する熱量は【大きくなる】
熱伝導率が小さいと移動する熱量は【小さくなる】
フーリエの法則
q  
dT
dx
熱流束
q
（※大きいか小さいかを問うているので、
スカラー量すなわち絶対値で考えよ！）
温度勾配
dT
[ K/m ]
dx
dQ
[W m 2 ]
dA
熱伝導率の影響
復習
q  
dT
dx
+
q
T2
T2
q
T1
より，熱伝導率λの単位は W/(m・K) となる
記入せよ
+
T1
熱伝導のしやすさ（熱の流れやすさ）の指標
7
復習
熱の出入りと温度変化の関係
8
熱の出入りがある棒の温度変化を考えてみる
復習
• 温度変化の基本原理
ある時点の温度T1, T2, T3が与えられたとき、Δt秒後の温度T2’を求めよ
T2’=T2+ΔT2 なので，
・熱を与えられるとそれに比例して温度が上昇
・熱を奪われるとそれに比例して温度が下降
ΔT2が分かればよい
この棒には，左からq12の熱量が
入り，右からq23の熱量が出ていく
熱量収支
Q  mc T
もらった
熱量
質量
T1
T2
T3
温度
上昇
比熱
Q  Q12  Q23
この熱量収支により棒の中点の
温度がΔT2変化する
Q  mcT2
・すべてはこの原理に基づいている
・「温度」という物理概念と、「熱量」という物理概念を
結び付ける法則＝原理
Q12
Δx
Q23
 mcT2  Q12  Q23
9
復習熱の出入りをフーリエの法則を使って表す
q  
dT
dx
①
①式を変化量で表すと
T2
T3
W = J/s
dQ
Q
T
q

 
dA tA
x
m2
④式と⑤式を③式に代入
 T T T T 
mcT2  Q12  Q23  tA 2 1  3 2 
x 
 x
  T  2T2  T1 
 T  2T2  T1 
 tA 3
  tA 3
 ⑥

x
x




s
T
x

T2  T1 
Q12    tA
④
x
T  T 
Q23  tA 3 2 ⑤
11
x
Q23
温度変化を求める
復習
Δz
 Q   tA
Δx
Q12
J
W/m 2
10
棒を三次元的に考えると
A
Δy
T1
③
A  yz
⑦
棒の質量は【密度ρ】×【体積】
Δx
m  xyz
⑧
12
2
温度変化を求める
復習
一次元熱伝導方程式の導出
復習
⑦式と⑧式を⑥式に代入
xyz cT2  t yz  T3  2T2  T1 
T1
Q12
T2
Δx
T3
x


 T3  2T2  T1 
 xcT2  t 

x


Q23
 T  2T2  T1 

 cT2  t  3
2
 x 

棒の中点の温度変化
を知ることができた！
T2 
t  T3  2T2  T1 


c  x 2

⑨
⑨式を変形すると
T2   T3  2T2  T1 

 ⑩

t
c  x 2

温度は一般的に位置と lim
時間の関数なので，
t  0
lim 温度の二階偏微分項を表す
x  0
 2T2 T3  2T2  T1

x 2
x 2
T T2

t
t
t  T3  2T2  T1 

 ⑨
 T2 
c  x 2

14
13
復習
一次元熱伝導方程式の導出
結果、⑩式は以下の様な偏微分方程式となる
離散式（温度変化）
T2
  T3  2T2  T1 



t
c  x 2

一次元熱伝導方程式
2
T
 T

t c x 2
式を比較してみよう
復習
⑪
極限化
lim
t  0 , x  0
温度Tの時間tによる微分
→ 温度の時間変化
比例する
つまり、温度勾配が大きく変化する場所では
温度の時間変化が大きい… ということ！
復習
離散化
（隣点の温度と分割幅Δx，
時間刻みΔtを利用）
T
  2T

t c x 2
温度Tの空間xによる2階微分
→ 「温度の勾配」の勾配
⑨
T1
連続式（熱伝導方程式）
Q12
15
数値シミュレーションの実際
温度T
時間t の変化によって…
方程式に現れる物理量は空間
および時間に対して 連 続 的
変換（離散化）が必要
5℃（固定）
17
Δx
Q23
16
実現象では
温度が連続的に変化
連続する時間を
Δtで離散化．
t0
t  t
t  2t 時間ステップ
t  3t
時間
コンピュータ上では物理量の 離 散 的
（デジタル的）な値しか扱えない
T3
時間の離散化
復習
• 対象となる現象を定式化
– 常微分方程式
– 偏微分方程式
T2
⑪
位置
5℃（固定）
シミュレーションでは，
時間的に離散（とびとび）
18
3
数値シミュレーションの限界
空間の離散化
復習
時間t=0 のとき
空間をΔxで離散化．
微小区間の温度を
代表点の温度で近似．
温度T
• とびとびの値を使って変化を模擬しても、極限的には
連続的な空間変化・時間変化へ近づけることが出来る
– 空間間隔を無限小にする・・・・・
– 時間分割を無限小にする・・・・・
lim
x 0
lim
t 0
位置
5℃（固定）
5℃（固定）
x
実現象では
温度が連続的に分布
節点
シミュレーションでは，
空間的に離散（とびとび）
・ 実際には限度がある
・コンピュータの計算能力
・ディスク容量
・メモリ容量
・計算時間の制限
19
20
偏微分方程式を離散的な式へ変換
変換対象は微分項
• 熱伝導方程式
T
  2T

t c x 2
• コンピュータで扱えるように，「連続的な」式（偏微分法
方程式）を「離散的な式」（差分方程式）に変換する
– 微小変量 → 有限の（具体的な）小さな値の変量
– 微分演算 → 割り算 に変換
– 積分演算 → 掛け算 ＆ 足し算 に変換
時間tに関する
一階偏微分項
• 変換は「近似」であるから，連続的な式との「誤差」が
生ずる
・・・ 一次元熱伝導方程式
空間xに関する
二階偏微分項
コンピュータでは扱えない
n 1
偏微分項を離散的に表したい
• 離散的な式を使うと，プログラム上の四則代数演算で
解くことが出来る
n
t
n 1
t
t
時間的・空間的な隣接点で表現できないか？
21
22
重要！
テイラー展開
• 偏微分の差分化は「テイラー展開」を用いて行う
• テイラー展開 (Taylor expansion) とは・・・
– ある関数の従属変数の微小変量に対する関数の値の変化
を、その関数の微係数の級数で表すこと
偏微分方程式の差分化
1   2u 
1   3u 
 u 
u x  dx   u x     dx   2  dx 2   3 dx 3  
2!  x 
3!  x 
 x 
1次の項
2次の項
3次の項
• 微小変量 dx → 有限変量 x に置き換えると
u j 1  u j  x
u
1
 2u
1
 3u
 ( x) 2 2  (x) 3 3  
x j 2!
x j 3!
x j
24
4
テイラー展開
テイラー展開の直観的な意味
• テイラー展開は、高次の項を考慮するほど漸近していく
• テイラー展開の直観的な意味
u j 1  u j  x
粗い近似 u j 1  u j
近似
1次の項
u j 1  u j  x
u j 1  u j  x
u
x
2次の項
変量 u
u：連続分布
u
x
j
u
1
 2u
1
 3u
 ( x) 2 2  (x) 3 3
x j 2!
x j 3!
x
・・
・
3次の項
u j 1
j
u
1
 2u
 ( x ) 2 2
x j 2!
x
u j 1  u j  x
高精度の
近似
u
1
 2u
1
 3u
 (x) 2 2  (x) 3 3  
x j 2!
x j 3!
x j
uj
x
x
j
節点番号・・・
j 1
高次の項ほど小さな値を取る
j
j 1
25
26
偏微分の差分化
一階の差分式の作り方
隣接する点で偏微分項を近似する
テーラー展開より、
• テイラー展開
u
1
 2u
1
 3u
u j 1  u j  x
 (x) 2 2  (x) 3 3  
x j 2!
x j 3!
x j
1次の項
2次の項
（その１）
3次の項
u j 1  u j  x
u
1
 2u
1
 3u
 (x) 2 2  (x)3 3  
x j 2!
x j 3!
x j
2次以降の項を無視
• 熱伝導方程式
T
  2T

t c x 2
時間tに関する
一階偏微分項
u j 1  u j  x
テイラー展開を用いて、偏微分項を
離散的な値で近似したい！
空間xに関する
二階偏微分項

u j 1  u j
j
x
u
 O( x 2 )
x j
u j 1  u j

 O( x)
x
u j 1  u j  x
u
x
j
隣り合う離散的な点の
物理量によって，一階の
微分項を近似
前進
j
j 1
後退差分 (Backward difference)
u
1
 2u
 ( x ) 2 2  
x j 2!
x j
u j  u j 1

 O ( x )
x
x
j 1
x
j
中心
j
j 1
– 前進差分と後退差分の差を取る
x
j 1
x
x
u
1
 2u
 (x) 2 2  
x j 2!
x j
差を取る
u
 2u
1
2
負方向 u j 1  u j  x x  2! (x) x 2  
j
j
Δxのかわりに
マイナスΔxを
用いれば…
2次の項まで考慮！
（但し相殺してゼロ）
u j 1  u j  x
u
x
（その3）
中心差分 (Central difference)
正方向 u j 1  u j  x
x
T
  2T

t c x282
近似によって
生じる誤差項
一階の差分式の作り方
（その2）
後退
一階の微分項を
近似できた！
 O( x)
27
前進差分 (Forward difference)
(First-order accurate
difference scheme)
隣接する点
u
x
一階の差分式の作り方
１次精度の差分法
u
 O( x 2 )
x j
3次以降の項を無視
2次精度の差分法
u
x
29

j
u j 1  u j 1
2 x
2
(Second-order accurate
difference scheme)
 O ( x )
30
5
一階の差分式の作り方
二階の差分式の作り方
（まとめ）
前進
後退
• 3種類の一階差分法
前進差分
u
x

• 二階の差分式
u  u (x )
– 一階の前進差分と後退差分の和を取る
1次精度
u j 1  u j
x
j
j
j 1
u
x
j 1
x
x
 2u
1
2
正方向 u j 1  u j  x x  2! (x ) x 2  
j
j
 O (x)
中心差分
和を取る
後退差分
前進差分
1次精度
u j  u j 1
u

 O ( x )
x j
x
後退差分
中心差分
u

x j
2
 u
2次精度
u j 1  u j 1
2x
2
T
 T

t c x312
2
 O (x )
x

2
時間に関する差分式
u j 1  2u j  u j 1
( x ) 2
j
後退差分
u
t
n

u n  u n1
 O ( t )
t
t0
t  t
t  2 t
t  3t
un
n 1
位置 x
位置x
u  u(t )
温度
温度
n は現在の時間ステップ
n+1 は次の時間ステップ
n-1 は前の時間ステップ
Δt は微小時間間隔
時間
時間
u n1
n
n 1
時間 t
１次元熱伝導の差分方程式（その１）
• シミュレーションに用いる差分式
– 時間に関しては・・・前進差分
– 空間に関しては・・・二階差分
n 1
n
t
j 1
j
x
n 1
n
T T j  T j

t 近似
t
n 1
t
j 1
x
T
  2T

t c x 2
t
・・・・・・・・・・ 各項の差分化
x
n 1
n
T T j  T j

t
t
t
二階差分
n
x 2

T jn1  2T jn  T jn1
x 2
n

n
n
 T j 1  2T j  T j 1
c
x2
・・・・・・・・・・ 未知数の分離
時間ステップ （※ 次数 ではない）
節点番号（空間分割）
 t
T jn1  2T jn  T jn1 
T jn 1  T jn 
c x 2
n
 2T T j 1  2T j  T j 1

x 2 近似
x 2
 2T
・・・・・・・・・・ 微分方程式の差分化
T jn 1  T jn
n
34
１次元熱伝導の差分方程式（その2）
• 熱伝導方程式を差分化する場合
前進差分
時間t
33
t
T
  2T

t c x 2
近似できた！
温度Tは場所xと時間tに
依存する関数 T(x,t)
温度T
時間ステップ （※ 次数 ではない）
u n1  u n

 O ( t )
t
 O ( x 2 )
時間的空間的に連続な温度分布
– 「時間的に隣接」する点との差分を考える
n
3次以降の項を無視
2次精度の差分法
T
  2T

t c x 2
• 時間の差分化
u
前進差分 t
 2u
1
x
j 1
j
j 1
u
2
負方向 u j 1  u j  x x  2! (x) x 2  
j
j
35
次の時間ステップn+1
（未知）
現時点 時間ステップn
（既知）
36
6
１次元熱伝導の差分方程式（まとめ）
• シミュレーションに用いる差分式
T jn 1  T jn 
次の時間ステップn+1
（未知）
 t
T jn1  2T jn  T jn1 
c x 2
現時点 時間ステップn
（既知）
T
  2T

t c x 2
コンピュータで扱える式になった！
現時点での値がわかれば
未来の値を計算できる！
n 1
コンピュータで
扱えない…
n
n 1
t
t
時間はΔtで離散化
次の時間ステップn+1に対して
前進差分
j
j 1
x
x
空間はΔxで離散化
両隣の節点j-1, j+1を用いて
中心差分
t
j 1
x
演習課題
37
演習1： Excelでの計算のやりかた
演習1： 1次元熱伝導
– 時間ステップ： Δt = 0.01
– 節点間の距離： Δx = 1
– c     1
20℃
壁の温度
（固定）
• 初期（t=0,n=0）において下図の様な温度分布で
あったとする．左右端の温度は変化しないとして
，真ん中の点の温度変化をExcelを用いて時間
ステップn=100 （時間 t=1.00） まで求めよ．
壁の温度
（固定）
5℃
5℃
j-1
j
20℃
壁の温度
（固定）
5℃
5℃
次の時間
ステップn+1
（未知）
j
j+1
下方向に時間軸、左方向に空
間軸をとったときの各位置の温
度Tjをまとめる
時間ステップ
j-1
 t
Tjn1  2T jn  T jn1 
c x 2
Δx
壁の温度
（固定）
 t
T jn1  T jn 
T jn1  2T jn  T jn1 
c x 2
T jn 1  T jn 
j+1
位置x の温度Tx [℃]
n
Tj-1
Δx
現在の時間ステップn
（既知）
Tj+1
5
20.00
5
1
5
?
5
2
5
?
5
3
5
?
5
…
…
…
時間
39
Tj
0
位置
計算対象であるTjは、過去時
点の温度Tj-1, Tj, Tj+1を用いて
算出する
40
演習1： Excelシートの例
演習1： データのまとめかた
重要！
1次元熱伝導の解析 ←1番上の行にタイトルを入れる習慣をつけよう（Sheet名にも）
 t
T jn1  T jn 
T jn1  2T jn  T jn1 数式も載せておけば、見直しが容易
（必須ではない）
c x 2
• 見やすくまとめること

– 他の人が見ても理解できるように、重要な情報を記載する
– 誰がいつ作ったのか？ 何のためのデータか？
– 計算条件はどれか？ 解（結果）はどれか？

解析条件
支配方程式
前進差分（時間）
中心差分（空間）
差分法
• 後で容易に変更できるように作ること
– 将来、計算式や条件を変更することを事前に想定し、
容易に変更できるよう工夫して作成する
– Excelでは絶対参照／相対参照などを有効に活用する
熱伝導方程式
41
位置x の温度Tx [℃]
Tj-1
Tj
Tj+1
5
20.00
5
1
0.01
5
19.85
?
5
2
0.02
5
19.70
?
5
3
0.03
5
19.55
?
5
4
0.04
5
19.41
?
5
5
0.05
5
19.26
?
5
6
0.06
5
19.12
?
5
7
0.07
5
18.98
?
5
温度固定
8
0.08
5
18.84
?
5
5
9
0.09
5
18.70
?
5
温度固定
10
0.10
5
…100まで
…
…
…
…
1
時間ステップΔt
0.01
比熱 c
1
密度 ρ
1
熱伝導率 λ
1
右端 j+1
t
0.00
境界条件
本講義・本課題のみの注意事項ではなく、今後のあらゆるデータ
整理・プログラムコード作成の際に気をつけておくべきこと
時間 [sec]
n
0
節点間距離 Δx
左端 j-1
時間ステップ
計算させるセルでは
5
18.57
?
解析条件セルを参照する
5
7
演習1： グラフも書いてみよ
時間的空間的に連続な温度分布
25
初期温度
時間
0.00
0.10
0.20
0.30
0.40
0.50
0.60
0.70
0.80
0.90
1.00
Temperature,℃
20
15
10
5
節点間の線は、
近似曲線で表示
0
0
0.5
1
1.5
温度T
温度Tは場所xと時間tに
依存する関数 T(x,t)
位置x
2
X
壁の温度
（固定）
壁の温度
（固定）
熱伝導の特性（温度の平滑化）が再現できる
時間t
44
一次元熱伝導方程式
演習2： 差分法
T   2T

t c x 2
• 正方向および負方向のテイラー展開より以下を求めよ
–
–
–
–
温度Tの空間xによる2階微分
→ 「温度の勾配」の勾配
温度Tの時間tによる微分
→ 温度の時間変化
45
ある場所における温度の時間的な変動量は，その場所に
おける2階微分値（とんがり具合・凹凸具合）に比例している
空間一階偏微分項に対する前進差分
空間一階偏微分項に対する後退差分
空間一階偏微分項に対する中心差分
空間二階偏微分項に対する二次差分
（講義中に説明した内容の復習）
• A4レポート用紙．手書き．来週授業開始時に提出
速い
遅い
46
47
本日のまとめ 【答案用紙に書いて提出】
• 微分方程式中の微係数を隣接する点で近似す
ることを ① 化という
②
• ①化には
展開を利用する
• 位置j点における物理量uのxに関する微分を以
下の３つの方式により表せ
u

③
– uj+1とujを用いて
x
u
– ujとuj-1を用いて x
– uj+1とuj-1を用いて u
x

j

j
④
⑤
j
• ③，④，⑤をそれぞれ，⑥ 差分，⑦ 差分，⑧ 差
分とよぶ
48
8