Transcript PowerPoint

熱伝導、波動のOpenGLによる可視化
11162791 日高 敬介
研究の概要、動機
• 概要
自然現象を記述する偏微分方程式
（熱伝導方程式,波動方程式）をいくつかの数値計算手法を
用いて計算する。またそれらをc言語の画像描画ライブラリ
OpenGLを用いて視覚化し、精度を比較する。
• 動機
福永先生からの提案,さらに自らも物理、微分方程式、アル
ゴリズムに関することを研究してみたいと思っていたから。
偏微分方程式
二次元熱伝導方程式
𝑢 𝑥, 𝑦, 𝑡 : 時刻𝑡での座標 𝑥, 𝑦 の温度、𝜶は定数
𝝏𝒖
𝝏𝟐 𝒖 𝝏𝟐 𝒖
=𝜶
+
𝝏𝒕
𝝏𝒙𝟐 𝝏𝒚𝟐
二次元波動方程式
𝑢 𝑥, 𝑦, 𝑡 : 時刻𝑡での座標 𝑥, 𝑦 の波の高さ、𝜶は定数
𝝏𝟐 𝒖
𝝏𝟐 𝒖 𝝏𝟐 𝒖
=𝜶
+ 𝟐
𝟐
𝟐
𝝏𝒕
𝝏𝒙
𝝏𝒚
上記の偏微分方程式を解くことで、熱の伝わりや水面の波の様子を表
現することができる。
しかし、方程式の厳密解を得るためには複雑な計算が必要であり、厳
密解も複雑になる。
偏微分方程式の離散化
偏微分方程式をその方程式に適した数値解析手法に基づき変形する。
（離散化）
代表的な数値解析手法
• 有限差分法
• 陽解法(Explicit_Method)
• 陰解法
• クランクニコルソン法(Crank_Nicolson_Method)
• 有限要素法(FEM)
有限差分法（陽解法）
差分法の陽解法について説明する。
一次元熱伝導方程式を陽解法で離散化してみる。
𝑢 𝑥, 𝑡 : 時刻𝑡での座標𝑥の温度
𝜕𝑢
𝜕2𝑢
=𝛼 2
𝜕𝑡
𝜕𝑥
微分の定義に基づいて、近似式に変形。
𝑢 𝑥, 𝑡 + ∆𝑡 − 𝑢 𝑥, 𝑡
右辺 ≅
∆𝑡
有限差分法（陽解法）
𝑢 𝑥, 𝑡 : 時刻𝑡での座標𝑥の温度
𝜕𝑢
𝜕2𝑢
=𝛼 2
𝜕𝑡
𝜕𝑥
𝜕𝑢 𝑥 + ∆𝑥, 𝑡
𝜕𝑢 𝑥, 𝑡
−
𝜕𝑥
𝜕𝑥
左辺 ≅ 𝛼
∆𝑥
𝑢 𝑥 + ∆𝑥, 𝑡 − 𝑢(𝑥, 𝑡) 𝑢 𝑥, 𝑡 − 𝑢(𝑥 − ∆𝑥, 𝑡)
−
∆𝑥
∆𝑥
=𝛼
∆𝑥
𝑢 𝑥 + ∆𝑥, 𝑡 − 2𝑢 𝑥, 𝑡 + 𝑢 𝑥 − ∆𝑥, 𝑡
=𝛼
∆𝑥 2
よって
𝑢 𝑥, 𝑡 + ∆𝑡 − 𝑢 𝑥, 𝑡
𝑢 𝑥 + ∆𝑥, 𝑡 − 2𝑢 𝑥, 𝑡 + 𝑢 𝑥 − ∆𝑥, 𝑡
=α
∆𝑡
∆𝑥 2
有限差分法（陽解法）
時間に関する漸化式の形に変形
𝑢 𝑥,𝑡+∆𝑡 −𝑢 𝑥,𝑡
∆𝑡
𝒖𝒕+𝟏
𝒙
=α
=
𝑢 𝑥+∆𝑥,𝑡 −2𝑢 𝑥,𝑡 +𝑢 𝑥−∆𝑥,𝑡
∆𝑥 2
𝒖𝒕𝒙
𝛼∆𝑡 𝒕
+ 2 𝒖𝒙+𝟏 − 2𝒖𝒕𝒙 + 𝒖𝒕𝒙−𝟏
∆𝑥
このような変形を行うことによりコンピュータに反復計算させるこ
とができる。 (上式では見やすさから 𝑢 𝑥, 𝑡 = 𝒖𝒕𝒙 とする)
離散式部分の実装例
有限差分法（クランクニコルソン法）
差分法のクランクニコルソン法について説明する。
𝜕2 𝑢
𝜕𝑥 2
を
𝜕2 𝑢(𝑥,𝑡+∆𝑡)
𝜕𝑥 2
と
𝜕2 𝑢(𝑥,𝑡)
𝜕𝑥 2
の平均で近似する。
途中計算を省き離散式を示すと、
−𝒖𝒕+𝟏
𝒙+𝟏
2∆𝑥 2
2∆𝑥 2
𝒕+𝟏
𝒕
𝒕+𝟏
+
+ 2 𝒖𝒙 − 𝒖𝒙−𝟏 = 𝒖𝒙+𝟏 +
− 2 𝒖𝒕𝒙 + 𝒖𝒕𝒙−𝟏
𝛼∆𝑡
𝛼∆𝑡
上式は等号より右側のすべての項が未知数であるので、陽解法
のように単純な漸化式にすることができない。
有限差分法（クランクニコルソン法）
𝒕+𝟏
そこで全要素(𝒖𝒕+𝟏
~𝒖
𝟎
𝒏−𝟏 )に対しての離散式を連立する。
−𝒖𝒕+𝟏
𝒙+𝟏
=
𝒖𝒕𝒙+𝟏
𝑎=
2∆𝑥 2
𝛼∆𝑡
+ 2, 𝑏𝑖 =
𝒖𝒕𝒊+𝟏
+
2∆𝑥 2
𝛼∆𝑡
− 2 𝒖𝒕𝒊 + 𝒖𝒕𝒊−𝟏 (1 ≤ 𝑖 ≤ 𝑛 − 1)
この連立方程式を解けば良いことになる。
2∆𝑥 2
+
+ 2 𝒖𝒕+𝟏
− 𝒖𝒕+𝟏
𝒙
𝒙−𝟏
𝛼∆𝑡
2∆𝑥 2
+
− 2 𝒖𝒕𝒙 + 𝒖𝒕𝒙−𝟏
𝛼∆𝑡
有限差分法（クランクニコルソン法）
連立方程式をコンピュータで解く手法はいくつかあるがここではLU分解
を用いた。
LU分解
係数行列𝐴を下三角行列𝐿と上三角行列𝑈に分解し、
効率よく掃き出し法（ガウスの消去法）を行う手法。
与えられた方程式
𝑨𝒙 = 𝑳𝑼𝒙 = 𝒃
に対し、変数 𝒚 を
𝑼𝒙 = 𝒚
とおき、これを上式に代入する。
𝑳𝒚 = 𝒃
から変数 𝒚 を求める。求めた解 𝒚 を 𝑼𝒙 = 𝒚 の右辺に代入し、解 𝒙
を求めることができる。
𝑳𝒚 = 𝒃はガウスの消去法の前進消去、𝑼𝒙 = 𝒚は後退代入に対応す
る。
有限差分法（クランクニコルソン法）
縦100*横200要素の二次元でのシミュレーションの実装を行ったため、
20,000元連立方程式となり、係数行列は20,000次正方行列となった。
最初に一度だけ行うLU分解は𝑛次正方行列に対して、計算量が𝑛3 に
比例し、アニメーションにするため繰り返しごとの計算量が𝑛2 に比例す
るので𝑛 = 20,000では計算に相当の時間を要する。
また、当開発環境では係数行列のための20,000* 20,000の配列のメモ
リの確保を行うことができなかった。
有限差分法（クランクニコルソン法）
今回現れる係数行列は
• 成分のほとんどが0である（疎行列）
• 0でない成分は対角成分付近に存在する（帯行列）
• 対称行列である。
以上の条件から確保するメモリを (帯の幅)× 𝑛 まで減らすことができた。
具体的には𝟐𝟎, 𝟎𝟎𝟎* 𝟐𝟎, 𝟎𝟎𝟎から𝟐𝟎, 𝟎𝟎𝟎* 𝟐𝟎𝟎となり、当初の1%の容
量で済ますことができた。
対称帯行列に対してのLU分解は「有限要素法概説 著者：菊池文雄」
を参考にした。
有限要素法
有限要素法について説明する。
方程式が定義された領域を小領域（要素）に分割し、各小領域におけ
る方程式を補間関数で近似する。
２次元でこの手法を用いる場合は主に三角要素を用いて分割する。
方程式をガラーキン法によって離散化する。ここで補間関数を 𝑁 とす
ると方程式にガラーキン法を適用した式は
𝑆
となる。
𝜕2𝑢 𝜕2𝑢
𝜕𝑢
𝑁 𝛼
+ 2 −
𝑑𝑆 = 0
2
𝜕𝑥
𝜕𝑦
𝜕𝑡
𝜕𝑢
𝜕2𝑢 𝜕2𝑢
=𝛼
+
𝜕𝑡
𝜕𝑥 2 𝜕𝑦 2
有限要素法
𝜕2𝑢 𝜕2𝑢
𝜕𝑢
𝑁 𝛼
+ 2 −
𝑑𝑆 = 0
2
𝜕𝑥
𝜕𝑦
𝜕𝑡
𝑆
左辺を 𝑘 と 𝑐 に分け、それぞれ計算すると、
𝑘 =
𝑆
=
𝑆𝑒
𝑐 =
𝑆
𝜕2𝑢 𝜕2𝑢
𝑁 𝛼
+ 2 𝑑𝑆
2
𝜕𝑥
𝜕𝑦
𝜕 𝑁 𝜕 𝑁 𝑇 𝜕 𝑁 𝜕 𝑁 𝑇
𝑁 𝛼
+
𝜕𝑥 𝜕𝑥
𝜕𝑦 𝜕𝑦
𝜕𝑢
𝑁
𝑑𝑆 =
𝜕𝑡
𝑆𝑒
𝑁 𝑁 𝑇 𝑑𝑆𝑒
𝑑𝑆𝑒
有限要素法(陰解法)
ここからは格子に対角線を引き三角要素分割したものを対象とすると、
−1 −1
𝛼 2
1
0
𝑘 =
2 𝑠𝑦𝑚
1
2
1 1
1
2 1
𝑐 =
24 𝑠𝑦𝑚
2
となる。これを元の式に代入すると、
𝜕𝒖
𝑘 𝒖 + 𝑐
=0
𝜕𝑡
離散化すると、
1
1
𝒏+𝟏
𝑘 +
𝑐 𝒖
=
𝑐 𝒖𝒏
∆𝑡
∆𝑡
有限要素法
三角要素分割では
２次元の有限差分法で
は、対象の点の上下左
右を含めた5点を材料
として、次ステップの
値を計算する。
一つ一つの三角要素に対して補間関数 𝑁 が作用するので接点5の値
を計算するために周りを含めた9点の値から次ステップの接点5の値を
計算する。
その為、一般的に差分法よりも精度が良くなる。
OpenGL
今回はプログラミング言語としてc++、
画像描画ライブラリとしてOpenGLを選択した。
OpenGL
非常に高速に動作し、画面に高精度な3D画像を描画できるライブラリ
関数例：
𝑔𝑙𝑉𝑒𝑟𝑡𝑒𝑥3𝑑𝑣(𝑐𝑜𝑛𝑠𝑡 𝐺𝐿𝑑𝑜𝑢𝑏𝑙𝑒 ∗ 𝑣)
→ 𝑣 3 = {𝑥座標, 𝑦座標, 𝑧座標} を引数とすると、指定した座標に点を
描画する。点の色、太さなどを設定できる。
他にも、自由に視点を設定できるなど、機能が充実している。
温度-色の関係式
熱伝導の可視化の際に温度を色合いで表現するためにサーモ
グラフィ風の関数を作成した。
温度によって色をなめらかに変化させるために色の変わり目部
分に三角関数を導入。
𝜋
sin(𝑥 − ) + 1
2
曲線部の関数 =
2
手法比較
• 初期条件
一辺100の正方形の上辺を100℃、その他を0℃。
𝒖(𝒊,𝟗𝟗) = 100 𝑖 = 0, … , 98 , 𝒖(𝟗𝟗,𝟗𝟗) = 0
𝒖(𝒊,𝒋) = 0 𝑖 = 0, … , 99, 𝑗 = 0, … , 98
𝝏𝒖
𝝏𝟐 𝒖 𝝏𝟐 𝒖
=𝜶
+ 𝟐
𝟐
𝝏𝒕
𝝏𝒙
𝝏𝒚
• 境界条件
一辺100の正方形の上辺は100℃、その他の辺は0℃。
𝒖(𝒊,𝟗𝟗) = 100 𝑖 = 0, … , 98 , 𝒖(𝒊,𝟎) = 0 𝑖 = 0, … , 99
𝒖(𝟗𝟗,𝒋) = 0 𝑗 = 0, … , 99 , 𝒖(𝟎,𝒋) = 0 𝑗 = 1, … , 98
• 終了条件
それぞれの手法で座標(50,50)の値のステップ間の差が
0.0001になった時。
𝒏
𝒖𝒏+𝟏
(𝟓𝟎,𝟓𝟎) − 𝒖(𝟓𝟎,𝟓𝟎) < 0.0001
• 比較対象値
𝒖(𝟓𝟎,𝒋) = 0 𝑗 = 5,10,15, … , 95
厳密解
手法比較
各手法の厳密解との差
0.25
FMEの厳密解との差
0.004
0.2
0.0035
0.003
0.15
0.0025
0.002
0.1
0.0015
0.001
0.05
0.0005
0
0
1
1
2
3
4
5
6
7
Explicit
8
9
10
11
12
Crank
13
14
15
16
FEM
17
18
19
2
3
4
5
6
7
8
9
10 11 12 13 14 15 16 17 18 19
手法比較
|old - new| < 0.0001 Explicit
Crank
FEM
Exact
Explicit
Crank
FEM
u[50][5]
1.71585
1.71584 1.75319 1.75304
0.03719
0.0372
0.00015
u[50][10]
3.47616
3.47614 3.54991 3.54961
0.07345
0.07347
0.0003
u[50][15]
5.32631
5.32629 5.43463 5.43419
0.10788
0.1079
0.00044
u[50][20]
7.31352
7.31348 7.45368 7.45312
0.1396
0.13964
0.00056
u[50][25]
9.48767
9.48763 9.65616 9.65552
0.16785
0.16789
0.00064
u[50][30]
11.9022
11.9021 12.0948 12.0941
0.1919
0.192
0.0007
u[50][35]
14.6147
14.6146 14.8265 14.8259
0.2112
0.2113
0.0006
u[50][40]
17.6873
17.6873 17.9132 17.9126
0.2253
0.2253
0.0006
u[50][45]
21.1871
21.1871 21.4213 21.4209
0.2338
0.2338
0.0004
u[50][50]
25.185
25.1849 25.4215 25.4215
0.2365
0.2366
0
u[50][55]
29.7542
29.7542 29.9873 29.9877
0.2335
0.2335
0.0004
u[50][60]
34.9673
34.9672 35.1909 35.192
0.2247
0.2248
0.0011
u[50][65]
40.8903
40.8903
41.099 41.1007
0.2104
0.2104
0.0017
u[50][70]
47.5755
47.5755 47.7639 47.7664
0.1909
0.1909
0.0025
u[50][75]
55.0501
55.0501 55.2136 55.2167
0.1666
0.1666
0.0031
u[50][80]
63.304
63.304
63.4384 63.4419
0.1379
0.1379
0.0035
u[50][85]
72.2772
72.2772 72.3791 72.3824
0.1052
0.1052
0.0033
u[50][90]
81.851
81.851
81.9179 81.9204
0.0694
0.0694
0.0025
u[50][95]
91.8479
91.8479 91.8781 91.8792
0.0313
0.0313
0.0011
Explicit Crank FEM
timestep
processing
time[sec]
12166 12168
32
153
3
1.7
手法の比較
有限差分法
手法
計算量
（繰返し時）
精度
*安定条件
任意形状
陽解法
クランク
ニコルソン法
𝑂( 𝑛 )
𝑂( 𝑛2 )
有限要素法
(陰解法)
𝑂( 𝑛2 )
低
条件安定
𝛼 < 0.25 (2𝐷)
適していない
高
無条件安定
適している
𝑢 𝑥, 𝑦, 𝑡 :
時刻𝑡での座標 𝑥, 𝑦
の温度、𝜶は定数
𝝏𝒖
𝝏𝟐 𝒖 𝝏𝟐 𝒖
=𝜶
+
𝝏𝒕
𝝏𝒙𝟐 𝝏𝒚𝟐
課題、展望
• 有限要素法の理解を深め、時間刻みと方程式を解く速度の関係に
関しての疑問を解消する。
• 有限要素法で任意形状のシミュレーションを行う。
• 境界要素値固定以外の境界条件でのシミュレーションを行う。
• より高度な事前現象に基づく偏微分方程式（KdV方程式など）の可
視化を行う。
KdV方程式
コルテヴェーグー
ド・
フリース方程式と
は浅いところでの
水の流れで生じる
波の現象を表した
方程式である。
参考文献
•
•
•
•
「数値計算工学」 著：森口繁一
「偏微分方程式の数値シミュレーション」 著：登坂宣好、大西和栄
「有限要素法概説」 著：菊池文雄
「流れと熱伝導の有限要素法入門」 著：矢川元基
• GLUTによる「手抜き」OpenGL入門
http://www.wakayama-u.ac.jp/~tokoi/opengl/libglut.html
• 明治大学数学科 桂田研究室HP
http://nalab.mind.meiji.ac.jp/~mk/program/
• Personal CAE Project
http://computation.cside.com/index.html