方程式解法ソフトによる
化学工学演習

「分離技術」 1999年1号 , pp.24-29に一部追加
新潟大学工学部化学システム工学科 伊東 章


ケミカルエンジニアの道具:計算尺からパソコンへ

化学プロセスシミュレータや3次元CADなど、現在のケミカルエンジニアの仕事はパソコン上の各種ソフトを自在に使いこなすことでおこなわれる。一方、彼らを育て送り出す役割の大学の化学工学の教育現場では、計算尺が電卓に進歩したものの、未だに図式解法や手計算、手書きの製図に終始している。これは化学工学の教科書自身がパソコンが道具となった現状に対応していないことが原因であるが、エンジニアの現場と教育との乖離がますます広がっているようである。

また、化学工学に限らず大学教育一般について最近は「魅力ある授業」が学生から求められている。化学工学演習では膨大な手計算を課すことになるが、もはや現在の学生に「計算力の修行」を強要することはむつかしく、化学工学科目の選択者を減らす一因となっている。

大学改革の流れのなか、化学工学の教育現場としても、現代のエンジニアの仕事法を反映し、かつ学生の手計算の負担を軽減した化学工学の授業・演習の新しいありかたを考えてゆく必要に迫られている。

BASIC or プロセスシミュレータ?

ひとむかし前はプログラミング言語「BASICによる化学工学演習」の構成が盛んに試みられた1,2)。しかし表計算ソフトの利用などで、エンジニアが自分でプログラミングをおこなうことはまれになった。これを反映してか、プログラミング言語による化学工学演習は普及しないまま止まっているようである。 

筆者もBASIC言語が搭載されたポケコンを学生に購入してもらい、BASICで化学工学演習をおこなっていたが、数年前に止めた。(その演習用に使っていたBASICプログラムは、当ページの「計算できる演習ページ」のVBScriptとしてしぶとく生き延びている。)

現在、ケミカルエンジニアの現場はプロセスシミュレータが大きな役割を演じている。

 

st991_55.gif (12553 バイト)

これは、パソコン画面上でプロセスを構成するだけで、物性値から収支計算までおこなってくれるパワフルな道具である。これを教育現場で使えれば教育効果がありそうにも思われた。しかし筆者の経験3)からすると、プロセスシミュレータはあまりにブラックボックス的で、これによる教育の効果には疑問を持っている。初級学生の目をひくためと、基礎は習得している修士課程学生の演習には有効だが、学部教育での化学工学演習に全面的に導入することは現実的ではないと考える。

ネットワーク上で利用する方程式解法ソフト

化学工学演習の要点は、基礎式である連立非線形方程式および微分方程式をたて、それを解くことである。図式解法は計算機の無い時代にこれらをなんとか解こうとした、先輩ケミカルエンジニアの偉大な工夫であった。しかし現在はパソコンを道具として方程式を直接解くことができる。

方程式解法ソフトとは、プログラミングをするのではなく、モデルを作成して解くべき式を書き下ろすだけでその方程式系を解析し、解を求めてくれるものである。プログラミング言語とプロセスシミュレータの中間に位置するものとして、方程式解法ソフトの化学工学演習への導入が教育効果の面で注目される。

なお、「Mathematica」のような「数式処理ソフト」と、ここで扱う「方程式解法ソフト」は類似しているが、「数式処理ソフト」は文字列操作(因数分解など)や数値計算からグラフまで、数学全般を扱うものであり、解析手続きを指示するための文法(コマンド)を覚える必要がある。これに対して「方程式解法ソフト」は「方程式を解く」ことに特化しており、式自身の記述はプログラミング言語風に式をそのまま書くだけで、解法の手順はソフトが自動解析してくれる。

方程式解法ソフトは以前から市販されていたものである4)が、価格のこともあり、電卓のように教室の学生全員が手元で使えるような環境はなかなか実現しなかった。やっと最近になり、化学工学分野生れで国産の方程式解法ソフトEQUATRAN-G5) (オメガシミュレーション)(以下EQG)にネットワーク対応版がでた。現在はどこの大学にもネットワークで結ばれたパソコン教室がある。そのような環境なら全端末でEQUATRANを使うことができる。

新潟大学の化学システム工学科ではLAN対応EQGを導入して、化学工学演習を試行中である。そのなかから例題を紹介する。

図2 学生各人が使えるLAN対応の方程式解法ソフト

方程式解法ソフト支援による化学工学演習

精留

図式解法の代表であるMcCabe-Thiele法は多変数の連立方程式を解く問題である。EQGによりこれを解くこと自身は簡単なので、実験との対比と操作変数を変えてのシミュレーションに演習の力点を置く。

演習  4段のオルダーショウ型蒸留器でメタノール/水系の精留を、還流比を変えておこなう。

物質収支から各段気液組成をあらわす連立方程式を書き、これを解いて実測の流出液組成xDと比較しなさい。

蒸気流量V [mol/s], 液流量L [mol/s]は一定と仮定する.全体の物質収支よりV =D+L, および zFV=DxD+Lx4が成立する.還流比(Reflux ratio)R をR=L/Dで定義する.8個の未知数 xD,, x1, y2, x2, y3, x3, y4, x4についての方程式は以下の8個である.

@全体の物質収支:
A塔頂から各段の下までの物質収支:
 
 
(式(2)〜(3)はyn+1xnの関係を与え,x-y 座標上では直線となる.これが操作線である.)
B各段を去る気液が平衡にあるという条件すなわち理論段の仮定より各段で次式が成立する.




(αは実測値(表1)を近似して  )

表1、図3にEQG上の記述と結果の図表示を示す。EQGの配列機能により多数の式をまとめて書くことができる。

表1 「精留」EQGソースリスト
VAR x(5), y(5),a(5)
zF=0.3; R=3
y(2:5)=(R/(R+1))*x(1:4)+(1/(R+1))*xD
a=4.609*x^2-9.886*x+7.748
y=a*x/(1+(a-1)*x)
y(1)=xD; y(5)=zF

図3 「精留」実験と計算の比較

吸収 

吸収も図解法を用いる代表的例題である。しかし,その図解法は積分ないし微分方程式の積分を図上でおこなっているものである。計算機の無い時代の工夫としてはすばらしいものであり、尊重されなくてはならないが、積分の道具があれば基礎式を直接扱うことも容易である。

演習  SO2を2vol%含んだ廃ガスを充填塔で水で洗浄することにより、出口濃度を1/20にしたい。 物質移動容量係数,kxa, kya [kmol/(m3 h)],ガス流量G [kmol/(m2 h)], 液流量L [kmol/(m2 h)]は表2のようである。この充填塔高さを設計しなさい。

吸収操作の代表的問題なので詳細は省略するが、普通はH.T.U. と移動単位数Nを導入して図式積分を利用する。ここでは基礎式の直接積分で解を求める。気相、液相のSO2モル分率をそれぞれy,xとした基礎式は以下のようである。

(1)

(2)

(3)

平衡関係: (4)

このy, xに関する連立常微分方程式を塔底(z=0)から積分する。初期値である塔底でのxの値,x1, は物質収支:
    G (y1-y2) =L (x1-x2)
からx1=0.000536 である。

以上をそのまま記述したのが表2、計算結果を図示したのが図4である。yiの入力を求められるが,yi=0.014など適当な値でよい。x の値が0になるz が解である。z =0.87 mとなった。

この方法は従来の教え方に比べて非常に簡単である。なお、H.T.U.やN が出てこないのは教育上問題があるかもしれない。

表2 「吸収」EQGソースリスト
y1=0.02; x1=0.000536
kya=309; kxa=5468; G=46.9; L=1659
G*y'=-kya*(y-yi)
L*x'=-kxa*(xi-x)
kya*(y-yi)=kxa*(xi-x)
xi=0.0329*yi+0.003735*SQRT(yi)
x # x1
y # y1
INTEGRAL z[0,2] step 0.01
trend x,xi,yi,y

図4 「吸収」計算結果 吸収塔内濃度分布

反応工学 

反応装置の解析も微分方程式が主役なので、方程式解法ソフトで演習ができれば演習そのものは「定式化」に集中でき効率的である。

演習6) 気相触媒反応の反応速度と反応速度定数:
  (5) (6)
が与えられている。触媒層入口からの触媒容積V に対して、転化率x と温度が連立常微分方程式:
(7)
(8)
で表せる。伝熱係数(UA/S)を変化させて計算をおこない、触媒層の温度が一定温度を越えない条件を探し、そのときの触媒層の容積Vを設計する。

式をそのまま記述するだけでよい。(表3)EQGのグラフ機能により計算結果をわかりやすく表示できる。(図5)

表3 「反応工学」ソースリストと計算結果
F0=40000; UAS=8E6
zA0=0.05; zB0=0.2; P=101.3; DA=(1-1-1)/1
Cp=31.4; Q=-142400; Tw=673; RB=1100
pA=zA0*(1-x)*P/(1+DA*zA0*x)
pB=(zB0-zA0*x)*P/(1+DA*zA0*x)
k=22400*exp(-75370/8.314/T)
r=-k*pA*pB
x'=-r*RB/(zA0*F0)
T'=(Q*r*RB-UAS*(T-Tw))/(Cp*F0)
x # 0
T # 673
INTEGRAL V[0,2] step 0.05
TREND x, T
output1 V,x,T step 0.05
---------------------------------------
計算結果例 (F0=40000; UAS=8e6)
---------------------------------------
V 1:x 2:T 
---------------------------------------
0 0000000 0.0000000 673.0000 
0.1000000 0.2320226 712.8071 
0.2000000 0.5260068 743.3337 
0.3000000 0.7411176 745.1796 
0.4000000 0.8430365 727.5061 
0.5000000 0.8902307 709.4480 
0.6000000 0.9160088 696.4883 
0.7000000 0.9325017 688.1395 
0.8000000 0.9442821 682.9569 
0.9000000 0.9533046 679.7632 
1.0000000 0.9605165 677.7752 

図5 「触媒層温度分布」計算結果の図示

混合 

トレーサー応答法による装置内混合過程の演習を実験装置と計算演習の組み合わせによりおこなう。

演習 完全混合、押し出し流れ、混合流れを模した装置内のトレーサーのインパルス入力にたいする出口での応答を測定する。

完全混合槽列モデル
(9)
および分散モデルによる理論応答曲線と比較して、混合のパラメータを決定せよ。

計算力の落ちている今の学生にとっては、(9)式のパラメータを変化させての繰り返し計算はたいへんな労力を要するようである。表計算なりEQGにより学生の計算労力への負担を軽くすることも、学習効率をあげることになるのではないか。レポート作成の過程でも計算に余計な労力を使わせないことにより、考察にも力点をおかせることとしたい。EQGによる記述(表5)と計算結果およびデータとの比較例を示す(図7)。EQGではREPEAT命令により繰り返し計算が容易である。

表4 「槽列モデル」のリストと計算結果
GLOBAL N=10
VAR x(N)
x(1)=1; x(2:N)=x(1:N-1)+1
F=PROD(x)/x(N)
E=(N^N/F)*t^(N-1)*exp(-N*t)
REPEAT t[0,2] step 0.25
TREND E

計算結果 N=10
------------------------------
t 1:E
------------------------------
0 0
0.250000 0.008629007
0.500000 0.3626558
0.750000 1.144405
1.000000 1.251100
1.250000 0.7651491
1.500000 0.3240717
1.750000 0.1065187
2.000000 0.02908153

図7 完全混合槽列モデルの理論応答曲線と実験との比較

制御 

プロセス制御や動特性の分野はラプラス変換など高等な数学が必要なため、教育上の工夫が必要である。ここで数学的処理を強調しすぎると学生の興味を遠ざけてしまいがちである。微分方程式を簡便に扱う道具があれば、複雑な数学上の手続きを経なくても、連立微分方程式の直接積分の演習によりプロセス制御・動特性の学習が容易になる。以下の実習例のように、積分の道具を前提とし、微分方程式だけで構成される新たなプロセス解析の教育課程も考えられて良いのではないか。

演習 この演習では液面制御実験と組み合わせておこなう。実験は、一定速度で流出しているタンクの液面をポンプで液を補給することで一定に保つもので、液面レベルを制御するためのポンプ流量を操作量とする比例-積分(PI)フィードバック制御系である1)。実験装置で限界感度法により制御パラメータを決定する。

次に、このプロセスが図のような制御系であるとして、数値モデル上で目標値rの単位ステップ状変化に対する制御量xおよび操作量uの応答をシミュレートしてみる。モデルにおけるの比例感度Kpと積分時間Tiを設定して制御の様子をみる。(下図の2.3は2.5に訂正)

モデルについて操作量u と制御量x の間の伝達関数を微分方程式で書くと次式。(’ は時間t 微分。)

(10)

調節計の動作も次の微分方程式で表せる。

 Ti u ' = Kp (Ti ε'+ε)  (11)

ここでとおくと、上式は、

  (12)。

t =0における初期条件は  x=x’=z =0、t= 0でr が0から1にステップ状に変化する。

EQGの記述と計算結果例を表5,図10に示す。EQGの自動グラフ表示機能により結果の図示も容易であり、パラメータを変えた試行が演習できる。さらに実験と対比することで理解が深まる。

表6 「制御シミュレーション」のソース
6*x''+7*x'+1=2.5*u
z=u-Kp*(r-x)
z'=Kp*(r-x)/Ti
Kp=1.8; Ti=1.5
x # 0
x' # 0
z # 0
r #1
INTEGRAL t[0,40] step 0.1
trend u,x
output1 t,u,x step 0.1

図 10制御モデルのグラフ表示

インターネット上の資料との連携

エンジニア教育の重要事項として、問題解決に必要な資料・データを自ら探すことの訓練がある。演習実施に際しては全てを問題文中に与えるのではなく、物性値などは便覧類から調査するように指定するのが望ましい。しかし、現実には全員に化学工学便覧を用意するわけにもいかない。そこで、学生実習用パソコンがネットワークに接続されていることを利用して、ネットワーク・インターネット上に化学工学便覧を置くことを試行している3)。(図11)このWeb上の便覧は全国どこからでも参照可能なので、どこの大学の演習にも利用できるし、化学工学教育の共通化にもつながる。

むすび

ここで紹介できなかった抽出、膜分離、カスケードなどの演習はWeb3)を参照されたい。例題をみておわかりのように、計算機上の演習こそ実験を同時におこなうなどして、数式と実際の現象との対応を把握させることが肝要である。

化学工学教育に危機意識を抱いている教官は過去、FORTRAN, BASIC, CAD, プロセスシミュレータなど最新の「道具」を教育に活用しようと試みてきている。しかしそれらも化学工学教育を根底から変革するには至らなかったようにみえる。ここで提言した方程式解法ソフト支援もまた結局はこれらの轍をふむことになるのかもしれないが、化学工学コースに優秀な学生をひきつけるために、現代的で魅力のある化学工学演習の開発に継続した努力をしてゆきたいものである。


参考文献

1) 化学工学協会編:”BASICによる化学工学プログラミング”, p.234, 培風館(1985)

2) 化学工学協会編:"化学工学プログラミング演習",培風館(1976)

3)伊東章:”化学工学資料のページ", http://irws.eng.niigata-u.ac.jp/~chem/itou/resource/res_home.html

4)伊東章:分離技術, vol. 17, p. 409 (1987)

5) http://www.omegasim.co.jp/product/eqg/indexj.htm

6)久保田宏、関沢恒男:"反応工学概論", p.149, 日刊工業新聞社 (1972)