PC 数学活用法
テキスト
静岡大学工学部物質工学科
木村元彦、前澤昭礼 著
1
第１章 VBA 活用法
目的
開発・製造現場で
大活躍
できる、全ての分野のエンジニアは、PC を用いた計算、計測、制御の技術を持っていなければならない。マイクロソフト社のエクセルには、非常に使用しやすい Visual Basic を基本としたプログラミング機能 VBA（Visual Basic for Applications）がある。 前期の 情報処理やデータ処理およびシミュレーションの授業で学習した VBA の機能をさらに活用できる 能力を修得することを目的とする。
1-1 ユーザーフォームを用いたプログラム
VBA を使用するために、「開発」→「Visual Basic」を選択し、VBE（Visual Basic Editor）の window を表示させる。
「挿入」→「ユーザーフォーム」を選択すると、図のようなユーザーフォームが作成される。また、同 時に、プロジェクトエクスプローラ内に、「UserForm1」が表示される。
図のように、ツールボックスから、CommandButton を２つ、TextBox を１つ、ユーザーフォームに作 成する。
2
CommandButton1 をクリックして、以下のプログラムを入力する。 Private Sub CommandButton1_Click()
TextBox1.Text = "こんにちは" End Sub
CommandButton2 をクリックして、以下のプログラムを入力する。 Private Sub CommandButton2_Click()
End End Sub これらのプログラムを入力した状態を図に示す。 フォームの表示とコード（プログラム）の表示を切り替えるためには、プロジェクトエクスプローラ内 の UserForm1 が選択されている状態で、 （コードの表示）あるいは、 （オブジェクトの表示） をクリックすればよい。
VBE の window の （実行）をクリックして、プログラムを実行してみる。CommandButton1 をクリッ
クすると、TextBox1 に「こんにちは」が表示される。エラーがあった場合には、 （リセット）をクリ
ックして、誤りを訂正する。
3
示）を選択すると、それぞれのシートに対応したコードの window が表示される。
すなわち、エクセルの VBA には、 各シート
This Work Book ユーザーフォーム
標準モジュール（後述する）
の４種類のプログラムを別々に記述できる場所があることを理解せよ。情報処理の授業のときに プログラムを記述したのは、Sheet1 に対応するコード window であった。
さて、この状態で、「開発」 → 「マクロ」としてみても、図のように、マクロは表示されない。つまり、 エクセル VBA では、ユーザーフォームは、エクセルの Sheet とは別の window であることが判る。
作成したユーザーフォームを、エクセルのファイルの読み込みと同時に表示できると便利である。
その場合には、プロジェクトエクスプローラ内の This Workbook を選択し、 （コードの表示）をク
リックする。そして、コードを表記する上の部分にあるオブジェクトボックス（左側）で Workbook を選 択し、プロシージャボックス（右側）で open を選択する。以下のコードを書く。
4 Private Sub Workbook_Open()
UserForm1.Show End Sub ファイルを保存し、エクセルを一旦終了してから、再度、ファイルを開いてみよ。作成したユーザー フォームが自動的に表示されることを確認せよ。 なお、ユーザーフォームを削除したい場合には、プロジェクトエクスプローラ内で、UserForm1 を右 クリックし、「UserForm1 の解放」を選択すると、作成した UserForm1 に関するフォームとコードの全 てを消去することができる。
1−2 標準モジュールを用いたプログラム
VBE の window を開き、「挿入」 → 「標準モジュール」を選択すると、図のような window が生成さ れる。
5 以下のコードを入力する。 Sub test() Cells(2, 3) = 123 Cells(2, 4) = "静岡大学" End Sub （実行）をクリックして、プログラムを実行し、エクセルの Sheet1 の中に”123”と“静岡大学”が 表示されることを確認せよ。 Sheet1 に表示されたものを消去してから、VBE のコードをつぎのように変更して、再度実行してみ よ。 Sub test() Sheet2.Cells(2, 3) = 123 Sheet3.Cells(2, 4) = "静岡大学" End Sub
今度は、Sheet2 と Sheet３の中のセルに表示される。このように、VBA では、全ての Sheet やユー ザーフォーム等に対して操作をすることができる。
Sub test() と End Sub までのプログラムを「プロシージャ」と呼ぶ。この状態で、「開発」 → 「マクロ」をクリックしてみる。すると、図のように、標準モジュール内に書かれたプロシージャは、エ クセルのマクロの一部になっていることが判る。1-1 で学習した、ユーザーフォームを使用したプロ シージャ以外は全てマクロとして表示される。（なお、Office2003 以前では、sheet に作成された VBA プログラムは、マクロとして表示されない。）
6
1-3 変数の扱い
ユーザーフォームを作成し、CommandButton を 2 つと TextBox6 つを設定し、以下のコードを入 力せよ。
Dim re
Private Sub CommandButton1_Click() vis = Val(TextBox1.Text)
7 dens = Val(TextBox2.Text)
vel = Val(TextBox3.Text) dia = Val(TextBox4.Text) re = dens * vel * dia / vis
TextBox5.Text = re End Sub
Private Sub CommandButton2_Click() Select Case re
Case Is > 4000: TextBox6.Text = "乱流" Case Is > 2300: TextBox6.Text = "遷移域" Case Else: TextBox6.Text = "層流" End Select
End Sub
（実行）をクリックして、プログラムを実行し、家庭用水道管に流れる水の Reynolds 数を求めて みる。（粘土：1mPa， 密度 1000kg/m3， 流速 0.1∼1m/s， 内径 0.013m）
正しく計算されることを確認したら、ファイルを適当なファイル名で保存してから、コードの先頭行 にある “dim re” を削除してから、再度実行してみよ。今度は、Reynolds 数は計算されるが、判 定が正しくなされないことが判る。この原因は、変数はプロシージャ内部でのみ有効なためである。 先頭のプロシージャの外部（上部）に、 dim re と定義をすることによって、re という変数は、全 てのプロシージャで共通に使用することのできる「グローバル変数」となる。プロシージャ内でのみ 有効な変数は、「ローカル変数」と呼ぶ。
8
定義してから実行してみよ。流速が 1m/s 未満の範囲では、Reynolds 数がゼロとなることを確認せ よ。
Private Sub CommandButton1_Click() Dim vel As Integer
vis = Val(TextBox1.Text) ・・・・・・・・ このように、変数に型を指定することができる。変数の型には、以下のものがある。数値計算な ど、計算精度が要求されるものには、Double（倍精度変数）を使用する必要がある。また、ループ 回数を指定する変数には、 整数型 Integer -32,768∼32,767 の範囲 整数値 Long -2,147,483,648∼2,147,483,647 の整数値 実数型 Single -3.402823E38∼-1.401298E-45 の負の値 1.401298E-45∼3.402823E38 の正の値 Double -1.79769313486232E308∼-4.94065645841247E-324 の負の値 4.94065645841247E-324∼1.79769313486232E308 の正の値 文字列型 String 0∼2GB の範囲のテキスト
例： Dim ABC As String * 25 25 文字以内の文字列 バイト型 Byte 0∼255 の整数値
ブール型 Boolean 真(True) または、 偽(False)
日付型 Date 西暦 100 年 1 月 1 日∼西暦 9999 年 12 月 31 日
通貨型 Currency -922,337,203,685,477.5808∼922,337,203,685,477.5807 の整数値 バリアント型 Variant データ型を明示しない（型が適宜変化する）
1-4 デバッグ
前項の Reynolds 数を計算する例題のファイルをロードし、UserForm1 に対するコードを表示させ よ。Private Sub CommandButton1_Click() 内にある dens = Val(TextBox2.Text) の行をクリック してから、ファンクションキー F9 を押す。これによって、この行にブレイクポイントが設定される。 UserForm1 を実行すると、ブレイクポイントで一時停止することを確認せよ。 動作が停止すると、コードの停止した行部分が黄色で表示される。このとき、コード内の変数 vis に マウスのカーソルを乗せてみると、vis に代入されている値が表示されることを確認せよ。 ファンクションキーF8 を押すたびに、一行ずつ実行され、実行された部分の変数の数値が変化す ることを確認せよ（ステップインと呼ぶ）。ブレイクポイントを解除するためには、ブレイクポイントの ある行内で、再度 F9 をクリックすればよい。 このように、ブレイクポイントとステップインを活用することによって、作成したコードの誤りを簡単 に見つけることができる。
9
1-5 DoEvents ・・・プログラムを暴走させない方法
エクセルをデザインモードとし、図のように、Sheet1 に CommandButton を２つと、Textbox を１つ 設定し、以下のコードを入力せよ。
Private Sub CommandButton1_Click() For i = 1 To 10000
DoEvents
Range("B2") = i Next i
End Sub
Private Sub CommandButton2_Click() End
10
デザインモードを終了して、CommandButton1 をクリックすると、Sheet1 のセルの数値が変化し、 10 秒間程度で停止する。再度、CommandButton1 をクリックし、その直後に CommandButton2 をク
11 リックすると、プログラムが終了することを確認せよ。 つぎに、コード内の DoEvents の直前（左側）に、「’」（半角のクォーテション）を付けて、コメント に変更せよ。（ コメント部分は実行されない。） 再度、CommandButton1 をクリックし、続 いて CommandButton2 をクリックしてみよ。今度は、停止させることができないことを確認せよ。 今、作成したプログラムは、For∼Next 文であるので、所定回数の後にプログラムは停止するが、 Do∼Loop 文を用いた場合には、無限ループとなり、システムを強制終了させる必要がある。 DoEvents は、繰り返しや sleep 文などの命令を使用する際に、プログラムを暴走から予防する有 効な方法である。また、エクセルでグラフを描画する場合にも、DoEvents を挿入することによって、 計算経過をグラフ表示させることができる。グラフ表示では、DoEvents を挿入していないと、最終 結果しか表示されない。
1-6 Function 文の使用法
Sheet1 に、CommandButton を１つ設定し、図のように、Sheet1 のセルにデータを入力せよ。次 に、Sheet1 に対応するコードウインドウに以下のコマンドを入力せよ。
Private Sub CommandButton1_Click() vis = Cells(3, 2)
den = Cells(4, 2) vel = Cells(5, 2) dia = Cells(6, 2)
Cells(7, 2) = re(vis, den, vel, dia) End Sub
Function re(a1, a2, a3, a4) re = a2 * a3 * a4 / a1 End Function
12 デザインモードを終了し、CommandButton1 をクリックすると、セル A7 に Reynolds 数が表示される。 このように、Function 文を使用すると、任意の関数を作成することができる。この例において、変 数 a1∼a4 を引数（ひきすう）と呼ぶ。また、Function からの計算結果を返値（かえりち）と呼ぶ。引 数や返値の変数型宣言をすることができる。（本来は宣言をすべきであるが、理解を簡単にする 目的でここでは宣言をしていない。）
1-7 マクロの記録の利用
エクセルには、「マクロの記録」という機能が付属している。複雑な操作に対するマクロを作成す る際、この機能は非常に便利である。図のように、Sheet1 に適当なデータを入力してから、「マクロ の記録」をクリックし、図のようにセルに色を付けたり枠を付けたりしてから、「記録終了」をクリック せよ。 セルの塗りつぶしと枠を元に戻してから、「開発」 → 「マクロ」と指定し、Macro1 を実行すると、 記録された内容が実行される。同様にして、Macro1 の内容を「編集」によって見ることができるば かりか、この標準モジュール内に自動的に作成されたプロシージャを適宜コピーして、他のプロシ ージャで利用することができる。この方法を用いれば、グラフ作成や印刷などエクセルの多くの複 雑な操作を、自作のプロシージャの中に使用することができる。また、マクロの命令を思い出せな いときにも、実際にエクセルのセルを使用して操作し、これを記録することによって、マクロの命令 を簡単に調べることもできる。14
レポート課題
以下に示すように、流速を ScrollBar で 0∼ 10m/s の間で変化させたとき、Reynolds 数を 自動計算して流れのパターン（層流、遷移域、 乱流）を判定するユーザーフォームとそのコー ドを作成せよ。 また、エクセルのファイルを読み込むと同時 にユーザーフォームが表示され、かつ計算を 連続実行するようにせよ。 レポート提出方法 レポート提出用の Web ページから提出するこ と。 http://ka31serv.eng.shizuoka.ac.jp/reportdata / ファイルは、Excel97-2003 形式（拡張子が xls） で保存し、ファイル名は、 学籍番号 8 ケタ ＋ “rep1” とせよ。例えば、50712001 の学生の場合には、 50712001rep1 となる。（拡張子まで表示され ていれば、 50712001rep1.xls である。） 提出期限 本日中（深夜 24 時まで）15
第２章 行列
目的 開発・製造現場で大活躍できる、全ての分野のエンジニアは、PC で行列を応用する能力が必 要である。会社での仕事（研究・開発）において、連立方程式の解法や、座標変換などでは、行列 を利用して計算をすると効率が非常に高くなることが多い。マイクロソフト社のエクセルには、行列 を計算する機能がある。この機能を知らないと損である。2-1 行列式の値
Sheet1 の適当な場所（図では、A2：C4 に 3 行 3 列の行列を作成せよ。（数値は適当で良い。） 別の適当なセルに、行列式の値を求める関数 =MDETERM(A2:C4) を入力せよ。行列の範囲指 定は、SUM などで範囲を指定するときと同様に、マウスで範囲選択をすることができる。2-2 転置行列
転置行列を求める操作は、少々複雑である。以下の操作を順番に実施すること。 ① セル A2 から C3 に適当な 2 行 3 列の行列を作成せよ。 ② セル E2 をクリックし、 =TRANSPOSE(A2:C3) と入力し、エンターを押す。すると、E2 のセルの 表示は、”#VALUE!”となる。このとき、行列の範囲はマウスで選択可能である。 ③ 転置行列が出力されるセル E2∼F4 を選択する。 ④ 数式バーをクリックする。（図示した状態） ⑤ 「Ctlr」キーと「Shift」キーを同時に押しながら、「Enter」キーを押す。 ← 重要！ この操作によって、行列が目的セルに展開される。16
なお、得られた行列を削除する際には、全てのセルを同時に削除する必要がある。（個別のセ ルは削除できない。）
17
2-3 逆行列
逆行列を求める操作は、転置行列を求める操作とほぼ同様である。（当然ながら、関数は異なる） 以下の操作を順に行う。 ① セル A2∼C3 に 3 行 3 列の行列を入力せよ。（ランダムな数値を入力せよ。） ② セル E2 をクリックし、 =MINVERSE(A2:C4) と入力し、エンターを押す。すると、E2 のセルに のみ逆行列の値が表示される。このとき、行列の範囲はマウスで選択可能である。 ③ 逆行列が出力されるセル E2∼G4 を選択する。 ④ 数式バーをクリックする。（図示した状態） ⑤ 「Ctlr」キーと「Shift」キーを同時に押しながら、「Enter」キーを押す。 ← 重要！ この操作によって、行列が目的セルに展開される。18
2-4 行列の積
行列の積は、 MMULT 関数を使用する。使用法は、TRANSPOSE や MINVERSE 関数と同様 である。図を参考にして、実施せよ。行列を展開する際には、同様に、「Ctlr」キーと「Shift」キーを 同時に押しながら、「Enter」キーを押す。
19
2-5 連立 1 次方程式を解く
逆行列と行列の積を使用すると、連立 1 次方程式の解を極めて簡単に求めることができる。こ こでは、例として、以下の連立 1 次方程式の解を求める。 8X + 0Y - 3Z = -1 2X - 1Y + 0Z = 0 -7X + 1Y + 2Z = 1 この連立 1 次方程式を、行列を用いて表現する。 8 0 -3 X -1 2 -1 0 Y = 0 -7 1 2 Z 1 これを AX＝B で表現すれば、解は、X＝A-1B となる。すなわち、連立一次方程式の係数の行 列の逆行列を求め、これに、連立一次方程式の右辺の行列との積を求めればよい。 図のように、連立一次方程式の係数と右辺の数値を入力し、解を出力させるセル（図では G2） に、 =MMULT(MINVERSE(A2:C4),E2:E4) を入力する。前項までと同様に、解の行列を展開せよ。 この例では、X=1，Y=2，Z=3 が解である。20
2-6 座標変換
行列の積を用いると、座標変換ができる。ここでは、一例として、原点を中心とした図形の回転を 試す。元の座標値（ｘ，ｙ）に対して、角度θだけ回転した座標（ｘ’，ｙ’）は、 ｘ’ cos(θ) ‒sin(θ) ｘ ｙ’ sin(θ) cos(θ) ｙ で求められる。 ① 図のように、Sheet の 2 行目と 3 行目に変換元となる座標値（ｘ，ｙ）の組み合わせを 10 組作 成せよ。図の例では、y = x2としている。セル A2 と A3 には、X と Y と項目を記入しておく。セ ル A5 と A6 も同様。 ② セル C10∼D11 に、cos(θ)あるいは、sin(θ)等に相当する式を入力せよ。 例えばセル C10 は、 =COS(B8/360*2*PI()) とせよ。 ③ セル B5 をクリックし、 =MMULT($C$10:$D$11,B2:B3) と入力し、Enter キーを押す。 ④ 数式バーをクリックしてから「Ctlr」キーと「Shift」キーを同時に押しながら、「Enter」キーを押し て、セル B5 と B6 に行列を展開する。 ⑤ セル B5 と B6 を選択し、その書式をセル C5∼L6 にコピーする。 ⑥ デザインモードで、Sheet1 に、２つの CommandButton を設定し、以下のコードを入力する。 ⑦ デザインモードを終了し、CommandButton1 をクリックし、図形が回転することを確認せよ。 ＝21 Private Sub CommandButton1_Click()
For i = 0 To 360 Step 10 DoEvents
Cells(8, 2) = i Next i End Sub
Private Sub CommandButton2_Click() End
22 レポート課題 最後の例題に改造を加え、適当な文字や図形を回転させるファイルを作成せよ。 レポート提出方法 レポート提出用の Web ページから提出すること。 http://ka31serv.eng.shizuoka.ac.jp/reportdata/ ファイルは、Excel97-2003 形式（拡張子が xls）で保存し、ファイル名は、 学籍番号 8 ケタ ＋ “rep2” とせよ。例えば、50712001 の学生の場合には、 50712001rep2 となる。（拡張子まで表示されていれば、 50712001rep2.xls である。） 提出期限 本日中（深夜 0 時まで）
23
第３章 フーリエ級数とフーリエ変換
目的 開発・製造現場で大活躍
できる、全ての分野のエンジニアは、信号処理や制御の能力が必要 である。交流信号を扱う際には、時間軸を基準に考えるだけでなく、周波数を基準に考えることの できる能力が必要となる。本章では、時間によって変化する信号を周波数成分として分析すること のできる、Fourier 級数および Fourier 変換をコンピュータ・シミュレーションによって視覚的に学習 することを目的とする。フーリエ変換は極めて基本的な概念であり、工学を学ぶ学生諸君は、フー リエ変換を理解していないと、就職してから会社での仕事で恥をかくことになる。 [１]Fourier 級数 時間と共に周期 2p で変動する関数 f(t)は、(1)式の Fourier 級数で表現することができる。これを、 実験的に確認してみる。 操作方法 ① エクセルのファイル”フーリエ級数.xls”を開く。マクロを含んでいるので、マクロを有効にする。 このプログラムは、(1)式のグラフを描画するだけの単純なものである。 ② (1)式の係数a
nおよびb
nとして、表１∼表 7 に示したパラメータを入力して、どのような周期波 形が合成できるか図示せよ。Reset でグラフ等を初期化できる。)
(
p
t
n
sin
b
p
t
sin
b
p
t
sin
b
p
t
n
cos
a
p
t
cos
a
p
t
cos
a
a
)
t
(
f
n n1
2
2
2
2 1 2 1 0·
·
·
·
+
+
·
·
·
·
+
+
+
·
·
·
·
+
+
·
·
·
·
+
+
+
=
p
p
p
p
p
p
表１ 波形合成パラメータ No.1 ｎ 0 1 2 3 4 5 6 7 8 9 10 cos 0 5 0 0 0 0 0 0 0 0 0 sin − 0 0 0 0 0 0 0 0 0 0 表 2 波形合成パラメータ No.2 ｎ 0 1 2 3 4 5 6 7 8 9 10 cos 0 0 0 0 0 0 0 0 0 0 0 sin − 0 5 0 0 0 0 0 0 0 0 表 3 波形合成パラメータ No.3 ｎ 0 1 2 3 4 5 6 7 8 9 10 cos 0 5 0 0 0 0 0 0 0 0 0 sin − 0 5 0 0 0 0 0 0 0 024 表 4 波形合成パラメータ No.4 ｎ 0 1 2 3 4 5 6 7 8 9 10 cos 3 0 -1 0 -0.22 0 -0.1 0 -0.06 0 0 sin − 2.5 0 0 0 0 0 0 0 0 0 表 5 波形合成パラメータ No.5 ｎ 0 1 2 3 4 5 6 7 8 9 10 cos 5 -2 0 -0.2 0 -0.09 0 -0.05 0 -0.03 0 sin − 0 0 0 0 0 0 0 0 0 0 表 6 波形合成パラメータ No.6 ｎ 0 1 2 3 4 5 6 7 8 9 10 cos 5 -3.173 0 1.013 0 -0.585 0 0.381 0 -0.256 0 sin − -0.313 0 0.313 0 -0.313 0 0.313 0 -0.313 0 表 7 波形合成パラメータ No.7 ｎ 0 1 2 3 4 5 6 7 8 9 10 cos 9.6 -0.313 -0.313 -0.313 -0.313 -0.313 -0.313 -0.313 -0.313 -0.313 -0.313 sin − -3.173 -1.571 -1.03 -0.754 -0.585 -0.468 -0.381 -0.312 -0.256 -0.209 [２]Fourier 変換 任意の周期関数 f(t)に対して、コンピュータを使用することによって、(1)式の係数
a
nおよびb
n を簡単に求めることができる。係数a
nおよびb
nは、(1)式から判るように、f(t)の周波数成分毎の 強度を表現したものである。すなわち、関数 f(t)は、周波数を変数とした係数a
nおよびb
nを用い て離散的に表現することができる。つまり、時間を変数とする周期関数 f(t)を、周波数を変数とす る関数 F(f)として表現することができることが推測できる。時間を変数とした関数 f(t)を、周波数を 変数とした関数 F(f)に変換することを Fourier 変換という。また、その反対に、周波数を変数とした 関数 F(f)を、時間を変数とした関数 f(t)に変換することを逆 Fourier 変換という。[１]で試した操作は、 逆 Fourier 変換を実施したことになる。 操作方法 ① エクセルのファイル”fft.xls”を開く。マクロを含んでいるので、マクロを有効にする。 ② 周期関数 f(t)をデータ番号０番から 31 番の 32 点のデータで離散的に入力する。入力する周 期関数として、表８に示したものを適用する。 ③ コマンドボタン“FFT”を押すと、FFT がマクロで実行されるので、計算結果を入力と共にグラフ で描け。25 表８ 関数 f (t) 1 cos(2πt) 2 cos(4πt) 3 cos(8πt) 4 sin(4πt) 5 cos(2πt) +sin(4πt) 6 sin(2πt) 0 < t < 0.5 秒 0 0.5 ≦ t < 1 秒 7 20t 0 < t < 1 秒 表９ 周波数 波形１ 波形２ 波形３ 波形４ 実数 虚数 実数 虚数 実数 虚数 実数 虚数 0 0 0 0 0 0 0 5 0 1 2 0 0 0 0 0 -2 0 2 0 0 2 0 0 2 0 0 3 0 0 0 0 0 0 -0.2 0 4 0 0 0 0 0 0 0 0 5 0 0 0 0 0 0 -0.09 0 6 0 0 0 0 0 0 0 0 7 0 0 0 0 0 0 -0.05 0 8 0 0 0 0 0 0 0 0 9 0 0 0 0 0 0 -0.03 0 10 0 0 0 0 0 0 0 0 11 0 0 0 0 0 0 0 0 12 0 0 0 0 0 0 0 0 13 0 0 0 0 0 0 0 0 14 0 0 0 0 0 0 0 0 15 0 0 0 0 0 0 0 0 16 0 0 0 0 0 0 0 0
26 [３]逆 Fourier 変換 周期関数 f(t)を構成する周波数成分（すなわち、周波数を変数とする関数 F(f)）を指定し、これか ら時間を変数とした関数 f(t)を求めることを、逆 Fourier 変換という。逆 Fourier 変換を幾つかの例 で体験してみる。 操作方法 ① エクセルのファイル”fft.xls”を開く。マクロを含んでいるので、マクロを有効にする。 ② 表９に示した周波数データを入力し、逆 FFT を実行する。コントロールボタン“逆 FFT”を押す と、逆 FFT がマクロで実行されるので、計算結果を画面にて確認せよ。 ③ 逆 FFT を実施した結果に対して、FFT を実施しても元の周波数データが変化しないことを確 認せよ。
FFT とは、 Fast Fourier Transform（高速フーリエ変換）の略であり、サンプル数が２のべき乗で あるときに Fourier 変換を極めて短時間に計算できるアルゴリズムである。サンプル数が２のべき 乗でない場合には、DFT；Discrete Fourier Transform（離散フーリエ変換）にて計算しなければなら
ない。データ数が N のとき、DFT ではおよそ N２回の計算が必要であるが、FFT では、（N/2）log ２N 回の計算ですむ。 レポート提出方法 全ての結果（グラフ）を適切な項目名と共に Word に張り付け、ファイルを提出する。提出先は、 授業用の Web 提出ページである。 http://ka31serv.eng.shizuoka.ac.jp/reportdata/ ファイル名は、 学籍番号 8 ケタ ＋ rep3 Word97-2003 の形式で提出すること。 提出期限 実習実施日を含めて２日以内 〔参考：Fourier 変換式の導出〕 (1)式の係数
a
nおよびb
nは、つぎのように求めることができる。 先ず、a
0 を求める。(1)式の両辺を t=0 から t=2p の範囲で積分する。)
2
(
sin
2
sin
sin
cos
2
cos
cos
2
)
(
2 0 2 0 2 2 0 1 2 0 2 0 2 2 0 1 2 0 0 2 0·
·
·
·
+
+
·
·
·
·
+
+
+
·
·
·
·
+
+
·
·
·
·
+
+
+
=
ò
ò
ò
ò
ò
ò
ò
ò
p n p p p n p p p pdt
p
t
n
b
dt
p
t
b
dt
p
t
b
dt
p
t
n
a
dt
p
t
a
dt
p
t
a
dt
a
dt
t
f
p
p
p
p
p
p
右辺は、その第１項を除き全てゼロとなるので、a
0 は、(3)式のようになる。27
)
3
(
)
(
1
2 0 0=
ò
pdt
t
f
p
a
つぎに、a
nを求める。(1)式の両辺にp
t
n
p
cos
を乗じてから、その両辺を t=0 から t=2p の範囲で 積分する。 ) 4 ( cos sin cos sin cos cos cos cos 2 cos ) ( 2 0 2 0 1 2 0 2 2 0 1 2 0 0 2 0 · · · · + + · · · · + + · · · · + + · · · · + + =ò
ò
ò
ò
ò
ò
p n p p n p p p dt p t n p t n b dt p t n p t b p t n a dt p t n p t a dt p t n a dt p t n t f p p p p p p p p p 右辺は、p
p
t
n
p=
ò
2 0 2cos
p
を含む項を除き全てゼロとなるので、a
n は(5)式のようになる。)
5
(
cos
)
(
1
2 0ò
=
p ndt
p
t
n
t
f
p
a
p
これは、n=0 のとき、(3)式のa
0と同値となるので、(
n
³
0
)
で有効となる。 つぎに、b
n を求める。(1)式の両辺にp
t
n
p
sin
を乗じてから、その両辺を t=0 から t=2p の範囲で 積分する ) 6 ( sin sin sin sin cos sin cos sin 2 sin ) ( 2 0 2 2 0 1 2 0 2 0 1 2 0 0 2 0 · · · · + + · · · · + + · · · · + + · · · · + + =ò
ò
ò
ò
ò
ò
p n p p n p p p dt p t n b dt p t n p t b dt p t n p t n a dt p t n p t a dt p t n a dt p t n t f p p p p p p p p p 右辺は、p
p
t
n
p=
ò
2 0 2cos
p
の項を除き全てゼロとなる。)
7
(
sin
)
(
1
2 0ò
=
p ndt
p
t
n
t
f
p
b
p
以上をまとめると、(1)式は、)
8
(
sin
cos
2
)
(
1 0å
¥ =÷÷
ø
ö
çç
è
æ
+
+
=
n n np
t
n
b
p
t
n
a
a
t
f
p
p
と書き換えられ、その係数a
nおよびb
nは、(5)式および(7)式で求められる。)
7
(
)
1
(
sin
)
(
1
)
5
(
)
0
(
cos
)
(
1
2 0 2 0³
=
³
=
ò
ò
n
dt
p
t
n
t
f
p
b
n
dt
p
t
n
t
f
p
a
p n p np
p
2
2
ix ix ix ixe
e
x
sin
,
e
e
x
cos
---=
+
=
を用いれば、これら(8)(5)および(7)式はつぎのように 表現できる。28
)
9
(
2
2
2
)
(
1 / / / / 0å
¥ =-÷÷
ø
ö
çç
è
æ
+
+
-+
=
n p t ni p t ni n p t ni p t ni ne
e
b
e
e
a
a
t
f
p p p p)
10
(
2
2
2
)
(
1 / / 0å
¥ = -÷
ø
ö
ç
è
æ
-
+
+
+
=
n p t ni n n p t ni n nib
e
a
ib
e
a
a
t
f
p p)
11
(
2
2
2
0 0 n n n n n nib
a
c
ib
a
c
a
c
=
=
-
-=
+
と定義すれば、)
12
(
)
(
p t ni n ne
C
t
f
på
¥ -¥ ==
)
13
(
)
(
2
1
dt
e
t
f
p
C
p t ni p p n p-ò
=
)
14
(
)
(
2
1
)
(
p t ni n p t ni p pf
t
e
dt
e
p
t
f
p på
¥ò
-¥ = --÷
÷
ø
ö
ç
ç
è
æ
=
(1)式から、n=1 の周期関数の周期は 2p（周波数は、1／2p）であることが判る。また、n=2 の周期関 数の周期は 2p／2（周波数は、2／2p）であることが判る。同様に、ｎ=n の周期関数の周期は 2p／ n（周波数は、n／2p）であることが判る。よって、p と n の変わりに、周波数ｆを用いて(6)式を表現す ることができる。すなわち、p
f
p
n
f
2
1
2
D
=
=
と書き換えると、(15)式となる。)
15
(
)
(
)
(
t
f
t
e
2dt
e
2f
f
i ft n ft i p p÷
ø
D
ö
ç
è
æ
=
å ò
¥ -¥ = -p p ここで、Df
®
0
とすれば、å
は、ò
で表現できる。また、p→∞となる。すなわち、)
16
(
)
(
)
(
t
f
t
e
2dt
e
2df
f
ò ò
¥ ipft i pft ¥ -¥ ¥ -÷
ø
ö
ç
è
æ
=
)
17
(
)
(
)
(
f
f
t
e
2dt
H
¥ -ipft ¥-ò
=
とおけば、)
18
(
)
(
)
(
t
H
f
e
2df
f
ò
¥ i pft ¥-=
すなわち、(17)式が Fourier 変換、(18)式が逆 Fourier 変換を表現する式となる。29
第４章 ＰＩＤ制御シミュレーション
目的 開発・製造現場で大活躍
できる、全ての分野のエンジニアは、装置を制御する能力が必要であ る。化学装置を含むほとんどの産業機器では、電子工学的な制御が使用されている。例えば、反 応器内の温度を一定に保ったり、流体の流量を一定に保ったりすることのできる能力が必要であ る。また、装置の立ち上げやシャットダウン時には、所定の温度変化や流量変化に従って温度や 流量等を制御しなければならない。ＰＩＤ制御は、このような制御を精度良く実現する方法である。 本章では、恒温槽内の油の温度をＰＩＤ制御するシミュレーションによって、制御法の概略を視覚 的に学習することを目的とする。 理論 機器の制御法として、大きく分類すると次の２種類に分かれる。 （１） open loop 制御法・・・・・・NC（数値）制御など（２） feed back 制御法・・・・・・ON-OFF 制御，PID 制御など
open loop 制御法は、制御したい変数（制御量という）（例：位置）を別の変数（例：時刻）によっ て、予め設定された関係式によって制御する方法である。例えば、パソコンのプリンターに使用さ れているステッピングモータは、期待どおりの角度だけ回転するものとして制御されており、もしも 期待した位置からずれていても（脱調しても）、そのずれが修復されることは無い。一方、温度制 御装置等において広範に使用されているｆeed back（負帰還）制御法では、制御したい変数（制御 量）（例：温度）を測定して、常にその目標値との偏差を小さくするようにコントローラの出力（操作 量）が制御されている。この例のように、feed back 制御法においては、制御量の目標値と実測値 との差を用いて操作量が常に修正される。 feed back 制御法には、いくつかの代表的な制御方法がある。 （１） ON−OFF 制御 （２） 比例制御（P 制御；Proportional control）
（３） 比例−積分制御（PI 制御；Proportional-Integral control）
（４） 比例−積分−微分制御（PID 制御；Proportional-Integral-Derivative control）
一例として、本章で扱うような、ヒーターによって温度を一定値に制御することを想定してこれら を説明すると、つぎのようになる。 ＯＮ−ＯＦＦ制御は、図１に示すように、制御量（温度）が目標値よりも小さい場合にはヒーター 電流をＯＮ、大きい場合にはＯＦＦとする、非常に単純な制御法である。具体的な例としては、バイ メタルを用いたサーモスタット（例：熱帯魚の水槽の水温制御、コタツの温度制御）がある。ＯＮ− ＯＦＦ制御では、制御量が目標値付近で振動し、その制御精度は、比較的低い。 比例制御は、目標値と現在値との偏差の大きさに比例した電力をヒーターに出力する方法であ る。比例制御を用いると、ＯＮ−ＯＦＦ制御よりも精度の高い制御ができる。しかしながら、図２に示 すように、比例制御においても、偏差をゼロにすることはできない。これは、ある程度の大きさを持
30 った偏差と操作量とがバランスしており、バランスを保つためには、ゼロでは無い偏差が必要であ るからである。また、比例制御の感度を極端に高く設定すると、ＯＮ−ＯＦＦ制御に近づき、制御量 （温度）が振動してしまう。比例制御における制御量の目標値 SV(t)、現在値 PV(t)、偏差 x(t)との 関係を(1)式のように定義する。
)
1
(
)
(
)
(
)
(
t
SV
t
PV
t
x
=
-
比例制御における操作量の出力（ヒーターへ通電する電力）y(t)は、(2)式で表現される。)
2
(
)
(
)
(
t
K
x
t
y
=
p pK
は、比例感度である。すなわち、制御量の偏差（目標温度からのずれ）に比例した操作量（ヒ ーター電力）を出力する方法である。なお、(2)式の右辺に、操作量のオフセット（初期値）を加える 方法もあり、適切なオフセットを与えると、比例制御における偏差を小さくすることができる。 ＰＩ制御は、比例制御の機能に加えて、制御量の偏差の時間的積分値を利用して、比例制御で は補正できない定常的に存在する偏差をほぼゼロにすることのできる制御方法である。ＰＩ制御法 を数式で表現すると、(3)式となる。ò
+
=
K
px
(
t
)
K
I tx
(
t
)
dt
(
)
)
t
(
y
03
IK
は、積分感度である。例えば、制御量（温度）が定常的に目標温度よりも小さい場合には、制 御量の偏差の積分値が増大し、操作量（ヒーター電流）を増大させることができる。 ＰＩＤ制御は、ＰＩ制御の機能に加えて、制御量の偏差の時間的微分値を利用して、外乱や目標 値の変化に対する応答性能を高めたものである。ＰＩＤ制御法を数式で表現すると(4)式となる。 DK
は、微分感度である。 ヒーターの例で、(4)式を言葉で表現すると、つぎのようになる。 擬人的な表現をすれば、それぞれの制御法は、つぎのように表現できる。 ＜Ｐ制御＞･･･現在の状態のみから判断して、現在値を補正する方法。 ＜ＰＩ制御＞･･･過去の反省と現在の状態から現在値を補正する方法。 ＜ＰＩＤ制御＞･･･過去の反省、現在の状態および未来の予測から現在値を補正する方法。 なお、積分感度K
Iは積分時間T
Iを用いて、 IT
Kp
と表現することもある。また、微分感度K
Dは微 分時間T
Dを用いて、KpT
Dと表現することがある。すなわち、一般的には、(4)式は、(5)式のよう に表現される。)
(
dt
)
t
(
dx
K
dt
)
t
(
x
K
)
t
(
x
K
)
t
(
y
D t I p4
0+
+
=
ò
(
比例感度)(
偏差)
＋(
積分感度)
(
偏差)
＋(
微分感度)
(
偏差)
＝ 出力パワー ｔ ０ dt d dtò
31 目標値 最大出力 時刻 温度 ﾋｰﾀｰ出力 目標値 最大出力 ﾋｰﾀｰ出力 温度 時刻 温度 ﾋｰﾀｰ出力 目標値 最大出力 時刻 ﾋｰﾀｰ出力 温度 目標値 最大出力 時刻 図１ ＯＮ−ＯＦＦ制御 図２ 比例制御 図３ ＰＩ制御 図４ ＰＩＤ制御
32
)
(
dt
)
t
(
dx
KpT
dt
)
t
(
x
T
Kp
)
t
(
x
K
)
t
(
y
D t I p5
0+
+
=
ò
操作方法 [１．定値制御] 図５に示すような、油の入った恒温槽の温度制御をシミュレーションする。演習では、(4)式に従 って、Kp
を比例感度、K
Iを積分感度、K
Dを微分感度として変化させてみる。 ① ファイル”PID.xls”を実行する。 ② 表１に示したパラメータにてシミュレーションを行い、その温度変化等を所定の用紙に記入せ よ。 表１ 操作パラメータ 比例感度 積分感度 微分感度 遅れ 初期温度 ゆらぎ 自動/手動 1 30 0 0 0 20 無し 自動 2 300 0 0 0 20 無し 自動 3 30 0 0 10 20 無し 自動 4 300 0 0 10 20 無し 自動 5 1000 0 0 10 20 無し 自動 6 30 50 0 10 20 無し 自動 7 30 50 0 10 90 無し 自動 8 − − − 10 20 ＯＮ 手動 9 30 50 0 10 90 ＯＮ 自動 10 30 50 200 10 90 ON 自動 下線（太字）で表現した部分が変更すべきパラメータである。恒温槽
PIDｺﾝﾄﾛｰﾗｰ
ﾋｰﾀｰ
電源
温度ｾﾝｻｰ
図５ 温度制御を想定する装置33 [２．プログラム制御] ① ファイル”PIDpro.xls”を実行する。遅れ時間が 10 秒になっていることを確認せよ。 ② 手動にて、画面に表示された温度変化曲線に従って温度を制御してみよ。 ③ 次に、自動（比例感度 50、積分感度 15、微分感度 50）にて制御してみよ。 レポート提出方法 全ての結果（グラフ）を適切な項目名と共に Word に張り付け、ファイルを提出する。提出先は、 授業用の Web 提出ページである。 http://ka31serv.eng.shizuoka.ac.jp/reportdata/ ファイル名は、 学籍番号 8 ケタ ＋ rep4 Word97-2003 の形式で提出すること。 提出期限 実習実施日を含めて２日以内 （参考：ＰＩＤ制御における伝達関数） (5)式は、時刻を変数として、制御量の偏差 x(t)と操作量 y(t)との関係を表現している。(5)式に対 してラプラス変換を適用すると、(5)式は周波数を変数として簡素に表現することができる。(5)式の 両辺を Laplace 変換する。 偏差 x(t)をラプラス変換すると、
X
(
s
)
制御変数 y(t)をラプラス変換すると、Y
(
s
)
となる。x(t) の積分項および微分項のラプラス変換に対しては、(6)および(7)式の公式を適用す る。)
(
)
(
x
)
s
(
sX
dt
)
t
(
dx
L
)
(
)
s
(
X
s
dt
)
t
(
x
L
t7
0
6
1
0-=
þ
ý
ü
î
í
ì
=
þ
ý
ü
î
í
ì
ò
すなわち、(5)式をラプラス変換すると(8)式となる。)
(
)
s
(
sX
KpT
)
s
(
X
s
T
Kp
)
s
(
KpX
)
s
(
Y
D I8
+
+
=
但し、x
(
0
)
=
0
とする 偏差X
(
s
)
と操作量Y
(
s
)
とは、伝達関数G
(
s
)
を用いて(9)式で表現できる。)
(
)
s
(
X
)
s
(
G
)
s
(
Y
=
9
(8)式を(9)式の形に整理すると、ＰＩＤ制御の伝達関数G
(s
)
が、(10)式となることが判る。)
(
s
T
s
T
Kp
)
s
(
X
)
s
(
Y
)
s
(
G
D I10
1
1
÷÷
ø
ö
çç
è
æ
+
+
=
=
34 ＜ラプラス変換＞ ラプラス変換対は、周期関数に限定される Fourier 変換対(11)をステップ関数などの非周期関 数に適用できるようにしたもので、(12)式で求められる。
df
e
)
f
(
H
)
t
(
f
ò
¥ i2pft ¥-=
H
(
f
)
¥f
(
t
)
e
-i2pftdt
(
11
)
¥-ò
=
ds
e
)
s
(
F
)
t
(
f
i st iò
-¥¥=
F
(
s
)
=
ò
¥f
(
t
)
e
-stdt
0s
=
i
2
p
f
(
12
)
ラプラス変換表
時間関数
ラプラス変換
)
(t
d
1
1
s
1
t
21
s
nt
1!
+ ns
n
ate
-a
s
+
1
(
)
t
ne
atn
-11
1
(
)
na
s
+
1
( )
w
t
sin
2 2w
w
+
s
( )
w
t
cos
2 2+
w
s
s
( )
t
e
-atsin
w
(
)
2w
2w
+
+ a
s
( )
t
e
-atcos
w
(
+
)
2+
w
2+
a
s
a
s
35
第５章 伝導伝熱シミュレーション
目的 開発・製造現場で大活躍
できる、全ての分野のエンジニアは、簡単な伝熱の計算ができなければならない。 特に、化学系と機械系エンジニアは、伝熱のシミュレーション技術も必要となる。多くの伝熱現象は非常に複 雑であり、コンピュータ・シミュレーションによって解析する方法が有効である。本章では、伝熱現象として最も 基本的である単純化した固体内の伝導伝熱をシミュレーションすることによって、伝熱現象の基本を理解す ると共に、コンピュータ・シミュレーションの基本的方法を理解することを目的とする。 理論・操作方法 ［１］ 定常１次元伝導伝熱 （先ず、１年次の復習から・・・・） 側面が断熱された、１次元の棒内を長さ方向に伝導伝熱する状態を考える（図１）。１次元の伝導伝熱を表 現する Fourier の式は、化学工学基礎で習得したとおり、(1)式で表現される。)
1
(
dx
dT
k
A
q
=
これは
、「単位時間、単位面積あたりの伝熱量(q/A)は、温度差（温度勾配）(dT/dx)に比例す
る。その比例係数は、熱伝導率 k である。」
ことを表現している。なお、マイナス符号は、xの増加する 方向に温度が増加する場合、伝熱方向が負の方向であることを意味していることも化学工学基礎で学んだ。 コンピュータで計算する為に、図１のように、棒を有限個数のセルに分割（離散化という）し、数値計算する ことを考える。分割された各セルの内部温度は、均一と考える。 図１ １次元伝導伝熱のモデル 側面が断熱された、１次元の定常伝熱では、温度分布は直線になることが直感的に理解できる。 （ここでは、厳密な導出は省略する。） 温度分布が直線になるということは、図１において、「注目セル T
iの温度は、左側セルの温度 T
i-1と
右側セルの温度 T
i+1の平均値になる」
ということを意味する。もしも、注目セルの温度が両側の平均値 になっていないと、左側からの伝熱量q
i-1と右側への伝熱量q
i+1に差が生じ、注目セルの温度 Ti が時間と 共に変化してしまう。（定常でなくなる。）T
iT
i-1T
i+1 側面断熱 1 + iq
1 -iq
注目セル36 Ti が、左右の温度の平均値になることを式で表現すれば、
)
2
(
2
1 1 + -+
=
i i iT
T
T
となる。 また、もしも、「端が断熱」とする必要がある場合には、「端部分の２つのセル温度が等しい」という条件設定 をすれば良い。すなわち、温度差が無ければ Fourier の(1)式から判るように、伝熱量がゼロとなる。 ＜操作方法＞ （先ずは、単純な例から開始） ① ファイル“定常１次元伝熱”を開く。シート１の説明を読む。 ② シート２を開く。計算方法を「手動」に変えよ。 （Office2003 以前では、ツール → オプション → 計算方法 → 手動）（Office2007 では、 Microsoft Office ボタン → Excel のオプション → 数式タブ → 手動）
③ 両端のセル以外のセル全てに、Ti=(Ti-1+Ti+1)/2 に相当する数式を代入せよ。 ④ 左端のセルに適当な温度（例えば 100℃）を入力し、右端のセルにも適当な温度（例えば 0℃）を入力せ よ。 ⑤ F9（再計算キー）を押し、棒内部の温度分布が収束するまで反復計算させ、収束した結果を図示せよ。 必要に応じ、再計算における計算回数の設定を変更すること。（Office2007 では変更が必要である。） ⑥ シート３を開き、同様に、両端のセル以外のセル全てに、Ti=(Ti-1+Ti+1)/2 に相当する数式を代入せよ。 ⑦ 右端のセルに、断熱の条件を入力せよ。すなわち、断熱の条件とは、「温度勾配がゼロである」ことであ る。具体的には、右端のセルとその隣のセルの値が常に等しくなるようにする。 ⑧ 左端のセルに適当な温度（例えば 100℃）を入力し、F9（再計算キー）によって収束値を求め、図示せよ。 ⑨ 終了時には、ファイルを保存しないこと。（保存する必要がある場合には、別名とすること。） ［２］ 定常２次元伝導伝熱 １次元の定常伝導伝熱は、直感的に理解できた。では、２次元では、どのようになるか考えてみる。２次元 の場合も、１次元のときと同様に、
「注目セル T
iの温度は、周囲のセル温度の平均値になる」。
２次 元の場合であっても、もしも、注目セルの温度が周囲の平均値でない場合には、注目セルへの伝熱量のバ ランスが崩れ、注目セルの温度が時間と共に変化してしまい、非定常となる。（図２） 式で表現すれば、)
3
(
4
1 , 1 , , 1 , 1 , + -+ -+
+
+
=
i j i j i j iJ j iT
T
T
T
T
となる必要のあることが直感的に理解できる。37 ＜操作方法＞ ① ファイル“定常２次元伝熱”を開き、説明書きを読んでから、sheet2 を開き、計算方法を「手動」に変えよ。 ② 厚さ方向に温度差の無い正方形の表面を持つ板状の固体内部の伝導伝熱を考える。板は、x，y それぞ れの方向に１１分割されて表現されている。 ③ 周囲を除く全てのセルに(3)式に相当する式を代入せよ。 ④ 周囲４辺のセルに一定温度（例えば０℃）を与えよ。 次に、板の中央付近の一点に一定温度（例えば 10０℃）を与えよ。 ⑤ Ｆ９キー（再計算キー）を用いて再計算を実行し、計算結果が収束するまで繰り返し、その結果を図示せ よ。（グラフが変化しないときは、グラフを選択し、右クリックし、「元データ」とし、列と行のチェックを入れ 替える。） ⑥ ④で一定温度（例えば 100℃）を与えたセルに０を代入してから、F9（再計算キー）を用いて再計算を実行 し全てのセルがゼロになるまで再計算させ、元の状態に戻す。 ⑦ 板の周囲の４辺に相当するセルに断熱条件を設定せよ。具体的には、図３に示すように、周辺のセルの 値がその中央側の隣りのセルと同じ値となるように設定する。角のセルには、その両隣りのセルの平均 値が計算されるように数式を入れよ。 ⑧ 板中央付近の１点（④で使用したセルと同一セル）に一定温度（例えば 100℃）を与えよ。 ⑨ Ｆ９キーを用いて再計算を実行し、計算結果が収束するまで繰り返し、その結果を図示せよ。 ⑩ 終了時には、ファイルを保存しないこと。（保存する必要がある場合には、別名とすること。） ⑪ なお、この計算において、グラフが動いて表示されるが、時刻による変化を示しているのではない。収束 値を求めるために、トライアンドエラーで反復計算しているので、グラフが変化していることに注意するこ と。
T
i,jT
i,j-1T
i,j+1T
i-1,jT
i,j+1 図２ ２次元定常伝導伝熱 定常ならば、Ti,jは、周囲の平均値38 セル C の温度は、セル A とセル B の平均値を代入せよ。 [３] 非定常１次元伝導伝熱 図１に示した棒内の非定常な伝導伝熱を考える。
D
t
秒間に、ｉ番目の注目セルの温度がD
T
iだけ上昇したと する。（図４参照）D
x
は、セルの間隔である。 図４ １次元非定常伝導伝熱のモデルt
D
秒間に、左側のセル（温度Ｔi-1）から注目セル（温度Ｔi）に移動した伝熱量は、(1)式よりt
x
T
T
kA
q
i i iD
×
D
-=
-)
(
1 1 (4) ・・・・（注目セルに入った量） ただし、Aは、伝熱面積、D
x
は、セル間の距離である。 一方、注目セル（温度Ｔi）から右側のセル（温度Ｔi+1）に移動した伝熱量は、 Ａ Ｂ Ｃ 図３ 断熱条件の設定T
iT
i-1T
i+1 側面断熱 1 + iq
1 -iq
x
D
At
x
T
T
kA
q
i i i×
D
D
-=
-)
(
1 1x
t
T
T
kA
q
i i i×
D
D
-=
+ +)
(
1 1t
T
xA
c
p×
D
×
D
i×
D
r
注目セル 入った量 たまった量 出た量39
