オメガシミュレーション

EQUATRAN-G

方程式解法ソフト

EQUATRAN-G（株）オメガシミュレーション

１．はじめに工学や物理などの分野に限らず、方程式によって問題を解くといったことは少なくありません。一般的に解析的な解が求められる方程式は少なく、数値計算によって方程式の解を求める方法が多く研究されてきました。技術計算の分野では技術者・研究者自らFORTRANあるいはBASICを使い、数値計算手法の知識を駆使してプログラムを作り、問題を解決してきました。しかし、満足するプログラムを作るまでには、多くの時間と労力が必要になっていました。

●数値計算の簡易言語
ここで紹介するEQUATRAN-Gはこのような技術計算の/*５元連立方程式*/A +B C + 2*D +E = 7
A + 6*B + 2*C - 3*D= 12
2*A B + 5*C D + 3*E = -2
-A - 2*B +C- 2*E = 3
A +B - 5*C +D - 4*E = 1
OUTPUT A, B, C, D, E
<<
A
B
C
D
E

計算結果

>>
= -1.264706= 3.738754= 3.624567= 5.472318= -2.794118

図1

５元連立方程式

悩みを一気に解決してくれるソフトです。技術計算をプロ
グラミングせずに簡便に行えるソフトで、「数値計算の簡易言語」と言えるものです。
EQUATRAN-Gは1985年に三井東圧化学が開発して販売をはじめたEQUATRAN-Mが原型になっています。
現在はWindows用としてEQUATRAN-Gを販売しています。
Windows XP/2000/NT/Me/98SE に対応しています。なお、大学向けの割引価格やネットワーク対応もありますので、お問い合わせください。

２．EQUATRAN-Gの概要●応用例
EQUATRAN-Gの基本となる機能は方程式の解法です。
方程式をそのまま入力すると直接数値解を得ることができます。たとえば、図1のような５元連立方程式を入力した
のち実行を指示すると、図の下のような計算結果を得ることができるのです。
しかし、単に方程式が解けると言われても、実際の仕事にどのように使えるのか、なかなか結びつかないかもしれませんので、そのあたりを説明しましょう。

1

オメガシミュレーション

EQUATRAN-G

一番身近な応用は、データから近似式を作るカーブフィ

します（デバッグします）。今度はその言語の文法などプッティングでしょう。観測データや文献などのデータから

ログラミングの知識が要求され、また多くの時間と労力が近似式を求めるのに最小2乗計算をし、データと近似式を必要となります。

グラフにプロットして評価できます。理系の学部を卒業しこれに比べ右側のEQUATRAN-Gの場合は、一番やっかた人なら必ずや経験があるでしょうし、実験を伴う研究業

いな部分であるアルゴリズムを考える段階と、プログラミ

務にはつきものといえるでしょう。

ングおよびデバッグの段階を自動化してくれます。数式モ

線形・非線形連立方程式を使えば、時間依存のない電気

デルを正しく入力しさえすればよいので、ほしいときに答

回路の計算、化学平衡計算、管路網の計算やプロセスの物えが得られると同時に、各人が抱えている本来の専門の仕質収支計算などのいわゆるバランス計算を解くことができ事に専念できることになります。

ます。
連続系のダイナミックシミュレーションは、常微分方程式でモデル化できます。振動系のシミュレーション、電気
回路の過渡現象の解析、制御系の解析や設計、反応速度の３．EQUATRAN-Gの各機能
つぎにEQUATRAN-Gの各機能をすこし細かく見ていきましょう。

検討、反応器のようなプラント機器の設計、血液循環系のシミュレーションや社会システムのシミュレーションなど●方程式解法機能
EQUATRAN-Gでは、線形連立方程式、非線形連立方程があります。

式、常微分方程式（高階または非線形を含む）、そして最●EQUATRAN とプログラミング言語との違い
プログラミング言語とはどこが違うのか、図 2 で説明しましょう。これは、問題を解く場合のアプローチの違いを示しています。
問

適化計算、最小2乗法計算（非線形を含む）を数値的に解くことができます。
方程式はそのままの形で入力すればよく、変形したり、解く順序に並び替える必要はありません。解く順序は、

題

問

数式モデルを作成する

EQUATRAN-Gが自動的に生成してくれます。

題

線形・非線形連立方程式の場合は、そのまま入力するだ

数式モデルを作成する

けです。特に非線形の場合は、一般に直接解くことができ計算のアルゴリズムを
考える

数値計算手
法の知識

プログラミングとデバ
ッグを行う

プログラミ
ングの知識

数式モデルを入力する

数値計算を行う

入力データ

数値計算を行うないので、繰り返し収束計算が必要になりますが、計算結果

FORTRANやBASICの場合

図2

EQUATRAN-Gでは線形部分や非線形部分を判断して、そ入力データ

れぞれの計算手法を自動的に組み込みます。なお、ユーザ

計算結果

ーが繰り返し収束計算の方法を指定することも可能です。

EQUATRANの場合

常微分方程式の場合は、積分の指定をする文をひとつ記述

問題解法のアプローチ

するだけで、高階の方程式でもそのまま扱うことができま

左側がFORTRANやBASICなどのプログラミング言語のす。最適化計算も独立変数と評価変数を指定するだけです。

場合です。数式モデルを作成したら、まず計算の手続き、

EQUATRAN-Gでは、以上の方程式を解く問題のほかに、

すなわちアルゴリズムを考えます。このとき、線形・非線

これらの方程式が混在した複合問題、たとえば常微分方程形方程式や常微分方程式などが含まれる場合は、それぞれ

式と線形・非線形連立方程式が混在した問題、最適化計算の数値計算の手法（たとえば、常微分方程式を解くのであの中に非線形連立方程式が含まれるような問題なども、扱

ればルンゲ・クッタ法など）の知識が必要になります。次うことができます。さらにユーザー関数を利用すれば、高にそのアルゴリズムにしたがってFORTRANやBASICによ度な複合問題（多重積分、2点境界値問題、動的システム

りプログラムを作成し、正しく計算できるように誤りを直のパラメータ同定問題、MINIMAX問題など）も扱えます。

2

オメガシミュレーション

EQUATRAN-G●方程式記述機能

表1

数式モデルを簡潔に記述するために、以下のような豊富な記述機能を持っています。
・配列変数‥‥1次元および2次元の配列変数が扱えます。

例題集［化学工学編］例題一覧

分類

タイトル

熱力学および物性

SRK状態方程式による圧縮係数の計算

熱力学および物性

Antoine式による蒸気圧の近似

熱力学および物性

定圧モル比熱の２次式への近似

熱力学および物性Wilson式のパラメータ推定および再現計算

熱力学および物性

NRTL式による液液平衡計算

熱力学および物性

化学平衡計算移動現象

層流境界層の速度分布の計算

粉粒体の特性

気流中の粒子の運動シミュレーション

流動

縮流管での圧力損失の計算

流動

管路網の流量と圧力の計算

流動コーナータップオリフィスの設計

流動

サージタンクの動的シミュレーション

伝熱

平板の１次元熱伝導の計算

伝熱

正四角柱の２次元熱伝導の計算

伝熱

向流型熱交換器の計算

伝熱

並流型熱交換器の計算

伝熱

1-2 shell-tube型熱交換器の計算

片対数・両対数グラフ、スプライン曲線による補間、1

伝熱

三連向流熱交換器の計算～3次式による近似曲線など、科学技術分野向きのグラフ

伝熱

三連熱交換器ネットワークの最適計算

蒸発三重効用缶の設計

蒸留

単蒸留の計算蒸留

平衡フラッシュ計算

・組み込み関数‥‥対数、指数、三角関数など36個の関数が内蔵されています。
・数表…変数間の関係が図や表として与えられている場合に、数表として定義し、方程式の中で関数のように利用できます。
・条件付きの式‥‥条件によって場合分けされるような式を表現できます。
・ユーザー関数とマクロ‥‥規模の大きな数式モデルはモジュール化して記述できます。

●グラフ作成機能機能が用意されています。自動設定機能により、簡単にグラフが作れ、しかも、各種のメイクアップ機能により、完成度の高いグラフが得られます。

●レポート作成機能
レポートの形に作成したテキスト（フォーム）中に、計算結果を任意の位置に、任意のフォーマットで埋め込むことができます。これにより、レポートの作成が容易にでき
ます。

蒸留

沸点計算

蒸留

簡易蒸留計算蒸留

McCabe-Thieleの図計算法による段数計算

蒸留

Thiele-Geddes法による蒸留計算

蒸留

回分蒸留計算蒸留

充填塔の充填高さと塔径の計算

吸収吸収塔の塔高計算

吸着・イオン交換

吸着平衡式のパラメータ推定

吸着・イオン交換

圧力スイング吸着の循環定常法計算

調湿・水冷却

空気調湿の計算

反応装置設計連続攪拌槽型反応器の解析

反応装置設計回分攪拌槽型反応器の解析

反応装置設計管型反応器（押し出し流れモデル）の解析

反応装置設計管型反応器（軸方向拡散モデル）の解析

おり、その応用分野は化学、機械、電気、制御などの各工反応装置設計管型反応器（半径方向拡散モデル）の解析

学分野はもちろんのこと、医学、薬学、経済学などにまで

生化学反応

酵素反応の追跡計算

プロセスの計画と設計

Williams-Ottoプロセスの物質収支計算

プロセスの計画と設計

Williams-Ottoプロセスの最適設計

プロセスの計画と設計

プロピレン製造プロセスの物質収支計算

プロセスの計画と設計

スチレン重合反応器のシミュレーション

プロセスの管理と制御

分解炉の温度制御シミュレーション

４．応用分野EQUATRANは現在国内で約2000セット以上が使われて及んでいます。数式モデルを用いる分野であればどんなところでも適用することができます。
特に化学工業分野では適応例が数多くあり、例題集にまとめられています。表1にその例題一覧を示しました。

3

オメガシミュレーション

EQUATRAN-G

５．例題
●例題１

二酸化炭素の熱容量（カーブフィッティング）

つぎのような二酸化炭素の熱容量の文献データがある。
温度をx、熱容量をyとして、最小２乗法で３次式に近似せよ。

y = a + bx + cx 2 + dx 3

温度[K]
250
300
400
500
700
1000
1200
1500

(1.1)

熱容量[J/mol･K]
35.44
37.52
41.44
44.68
49.59
54.32
56.34
58.38

図3

EQUATRAN-Gのウィンドウ

問題を解くために数式モデルを入力して、ソーステキストを作成します。新規に作成する場合は[新規作成]コマ
ンドを、すでにあるファイルを修正するのであれば、[開・数式モデルく]コマンドを実行します。子ウィンドウとしてエディタ(1.1)式の係数 a、b、c 、dを決定するわけですが、これらにより計算される熱容量の値をycalとして、i番目の測定が開かれますので、図4のようにここにソーステキストを入力していきます。

値に対する量に添え字iを付けて表すとycali = a + bx i + cx i + dx i
2

3

(1.2)

となります。最小２乗法の考え方により、誤差eの２乗和S

ei = y i − ycali
S=

(1.3)

8

∑e

2

(1.4)

i

i =1

を最小とする a、b、c、d を求めることになります。(1.1)
式は係数について線形になっていますので、正規方程式を導出して解けば一意的に決めることができます。ここでは別の方法、すなわち最小２乗法問題を解法する機能を使っ

図4

例題１のソーステキスト

2て解いてみることにします。なお、a のように係数にべき乗が付いていたり、sin a のように係数が初等関数のなかにあるなど非線形の場合には、正規方程式が利用できません
ので、こちらの方法を使うことになります。
・ソーステキストの作成と計算の実行
EQUATRAN-Gを起動してみると、図3のようなウィンド
ウが開かれます。

/*と*/で囲んだものはコメントであり、１行目は特別にタイトルとして扱われます。
3行目のGLOBAL文はNというパラメータの定義をしており、配列変数の配列の大きさを与えています。5～8行目のVAR（VARIABLEの略）文は変数の定義をする文です。配列変数は使用する前に定義する必要があります。一方、ス
カラー変数はBASICと同じように定義を省略することもで4

オメガシミュレーション

EQUATRAN-G

きますが、定義しておけば変数の後ろに付けた説明項が表

EQUATRAN-Gには自動設定機能があるので、縦軸・横軸

示の際に使われますし、後からこのソーステキストを見たに変数名を指定するだけで、とりあえずグラフができあがときに分かりやすくなります。

ります。あとは自分の気に入るように編集してやりますが、

10～13行目は配列変数 x と y に実験データの値を与え対応するところ（たとえば、軸の上下限値を変更するのでています。16～17行目は、(1.2)、(1.3)式がそのままの形であれば、軸の近傍）でマウスをダブルクリックすると、修記述されているので説明はいらないでしょう。

正用のダイアログボックスが表示されますので、簡単に編

そして、19行目のFIND文が最適化計算あるいは最小２乗法計算を指示する文で、この場合、LEASTのつぎにある e

集ができます。最終的に完成したグラフは図6です。この３次式は文献データをよく表現できていると言えます。

の２乗和を最小にする a、b、c 、dの値を求めることを指なお、EQUATRAN-Gのグラフ機能は、片対数・両対数

示しています。また、#1 は初期仮定値を、UNTIL項は計グラフ、スプライン曲線による補間や、直線、２次・３次

算精度を指定しています。21行目のOUTPUT文は計算結果

曲線近似などの科学技術グラフ向けの機能をもっています。

を出力する変数の指定です。
[ラン]コマンドを実行して計算の実行を指示すると、計算結果ウィンドウが開いて、図5のような結果が得られました。
したがって、近似式は

y = 22.23 + 0.0621x − 3.923 × 10−5 x 2+ 9.2670 × 10− 9 x 3(1.5)ということになります。

図6

●例題２

例題１の計算結果（グラフ）アンモニア製造プロセスの物質収支（静的なバランス計算）

アンモニアは、以下のように窒素と水素から直接合成される（Haber 法）。
N2 + 3H2
図5

例題１の計算結果

・グラフの作成

2NH3

設計の第１段階としてプロセスの大まかなバランスを計算することとし、図 7 のような簡略化したフローを考

つぎに求まった係数が妥当なものか、グラフにして評価える。原料は乾燥した水素と窒素ガスで、不純物としてしてみます。先ほど実行したときに計算結果が出力されたアルゴンとメタンを含む。反応器では 200atm、500℃の条

リザルトファイルを指定して、[実行]-[グラフ]コマンドを件で触媒の存在下で反応が行われる。反応器出のガスは実行します。グラフには縦軸・横軸の上下限値や、目盛や冷却、凝縮により分離し、未反応ガスはリサイクルさ

罫線の間隔など多くの設定する項目がありますが、て、液体アンモニアは製品として取り出される。また、5

オメガシミュレーション

EQUATRAN-Gが成り立ちます。なお、不純物のAr、CH4は反応に関与し未反応ガス中の不純物の蓄積をパージにより防ぐ。
原料の組成（モル分率）がおのおの N2 : 0.240、H2 :ませんので、入ったものがそのままの流量で出ていきます。

0.743、Ar : 0.006、CH4 : 0.011 で、その流量が 100 [kmol/h]f 2 , Ar = f 3, Ar(2.5)

f 2 ,CH 4 = f 3,CH 4(2.6)の場合の物質収支計算をせよ。
化学平衡定数 K を組成を xj、圧力を P として以下のよ
うに定義し

分離器まわりでは、
2x NH 3x N 2 xH 2 3

K = K ' P2 =f 3, j = f 4 , j + f 5, j

(2.1)

(2.7)

であり、パージする分配器のまわりでは、パージ比率がαK=0.35 の値を用いる。分離器での液体アンモニアへの分ですから

離は、入ってくる成分流量に対する分離比で与え 0.512とし、未反応ガスのパージ比率はα=0.02 とする。また、各ガスのアンモニアに対する溶解度[kmol/kmol]は、それ

f 5, j = f 6, j + f 7 , j(2.8)

f 7 , j = f 5, jα(2.9)

ぞ れ N2 : 0.00160 、 H2 : 0.00740 、 Ar : 0.00043 、 CH4 :が成り立ちます。

0.00225 とする。1)

そして、ストリームのトータル流量をFi [kmol/h]とする
未反応ガス

と
6

7
5

原料

1

2

反応器

3

パージガス

Fi =

5

∑f

i, j

(2.10)

j =1

分離器

で計算できます。
4

図7

アンモニア

アンモニア製造プロセス

以上の関係をソーステキストにしたのが図8です。３行
目のLOCAL文はパラメータを定義しており、実行時にこれらの左辺の文字列は右辺の数値に置き換えられます。10行

図7のように各ストリームに１から番号を振り、各スト

目は原料の成分流量です。そして、13行目が(2.2)式、16～

リームについて、N2 、H2 、NH3 、Ar、CH4 の５成分の流量

19行目が(2.3)～(2.6)式、24、32、33、35行目が(2.7)～(2.10)

を考えます。すなわち、ストリーム数(7)×成分数(5)の２次式であることは容易に分かると思います。

元の配列変数である、成分流量 fi,j [kmol/h]を定義して用い

20～21行目が化学平衡式(2.1)式ですが、組成 xj はfi,j / Fiます(i=1,2...7、j=1,2...5）。

で計算しています。また、この化学平衡式は非線形になっ

反応器前の混合器では、ており、解くために収束計算が必要になります。そのため

f 1, j + f 6, j = f 2 , j(2.2)の 繰 り 返 し 収 束 計 算 の 指 定 が 35 行 目 に し て あ り ま す 。
f(3,H2)を繰り返し変数に、収束のチェックをする式にeqのが成り立ちます。

ラ ベ ル を 付 け た 20 行 目 の 式 を 指 定 し て い ま す 。 な お 、

反応器まわりでは、化学量論からN原子およびH原子につf(3,H2)の後の数字はそれぞれ順に初期仮定値、下限値、上いて、入り出の原子数が等しいので

f 2 , N 2 × 2 + f 2 , NH 3 = f 3, N 2 × 2 + f 3, NH 3

限値を意味します。
(2.3)

f 2 , H 2 × 2 + f 2 , NH 3 × 3 = f 3, H 2 × 2 + f 3, NH 3 × 3(2.4)

25～29行目は、液体アンモニア中の各成分流量を、アンモニアの分離比と各ガスの溶解度を使って与えています。
39行目は出力する変数の指定ですが、f~はfの第１添字と6

オメガシミュレーション

EQUATRAN-G

第２添字を入れ替える転置配列を指定しています。これにここでは簡単な例を示しましたが、実務ではストリーム

より横方向にストリームを縦方向に成分を出力して見やすの数が50～100、成分数も20程度は扱うことになりますし、

くしています。

熱収支を合わせて計算することもあります。このような、図9に計算結果を示しました。トータル流量と成分流量が求まりました。
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:
27:
28:
29:
30:
31:
32:
33:
34:
35:
36:
37:
38:
39:

/*

大規模な非線形連立方程式を解く問題にも、EQUATRANGを適用することができます。

アンモニア製造プロセスの物質収支

*/

●例題３

LOCAL N2=1, H2=2, NH3=3, AR=4, CH4=5
VAR f(7,5)"成分流量[kmol/h]",F(7)"トータル流量 [kmol/h]",alpha = 0.02 "パージ比率[-]",K= 0.35 "化学平衡定数"

制御系の応答（ダイナミックシミュレーション）

図10のようなフィードバック制御系で、目標値 ys が単/* 原料 */
f(1) = ( 24.0, 74.3, 0, 0.6, 1.1 )位ステップ状に変化するとき、制御量 y および操作量 uの応答を求めよ。PIコントローラのゲインKと積分時間Tiは/* 混合器 */f(1) + f(6) = f(2)それぞれ/* 反応器 */f(2,N2)*2 + f(2,NH3)= f(3,N2)*2 + f(3,NH3)
f(2,H2)*2 + f(2,NH3)*3 = f(3,H2)*2 + f(3,NH3)*3
f(2,AR) = f(3,AR)f(2,CH4) = f(3,CH4)
eq: (f(3,NH3)/F(3))^2= K * (f(3,N2)/F(3)) * (f(3,H2)/F(3))^3

Ti = 5.2

K = 1.8,とする。また、操作量は弁開度であって|u|≦ 0.5 という制限があるとする。2)/* 分離器 */f(3) = f(4) + f(5)f(4,NH3) = f(3,NH3) * 0.512
f(4,N2) = f(4,NH3) * 0.00160
f(4,H2) = f(4,NH3) * 0.00740
f(4,AR) = f(4,NH3) * 0.00043
f(4,CH4) = f(4,NH3) * 0.00225

ys +

⎛
1 ⎞
⎟⎟
K ⎜⎜1 +
T
is ⎠
⎝

e

-

図10
/* 分配器 */f(5) = f(6) + f(7)f(7) = f(5) * alpha;2.5(1 + 3s)(1 + 2 s )(1 + s)y

フィードバック制御系

制御対象の伝達関数はF = SUM( f )
RESET

u

Y ( s)
2. 5=
U ( s) 6s 3 + 11s 2 + 6s + 1

f(3,H2) # 500[0,1000] BY eq

OUTPUT F, f~(3.1)

ですが、これを微分方程式の形に直すと

図8

例題２のソーステキスト

6
/*
<<
F

アンモニア製造プロセスの物質収支 */
計算結果
>>
=
: トータル流量 [kmol/h]1) 1002) 770.0656
3) 726.7461
5) 683.7404
6) 670.0656
7) 13.67481
f~=: 成分流量[kmol/h]
( 1)
( 2)
( 3)
1) 24
135.3393
113.6796
2) 74.3
515.6024
450.62313) 039.70623
83.02574
4) 0.629.10433
29.10433
5) 1.150.31336
50.31336
( 5)
( 6)
( 7)
1) 113.6116
111.3393
2.2722322) 450.3085
441.3024
9.0061713) 40.51656
39.70623
0.8103312
4) 29.08605
28.50433
0.5817211
5) 50.21772
49.21336
1.004354

図9

例題２の計算結果

d 3y

d2y

dt

dt 2

+ 11
3

+6

dy
+ y = 2. 5u
dt

(3.2)

となります。
4) 43.00568( 4)0.06801468
0.3145679
42.50918
0.01827895
0.09564565

また、PI制御器の伝達関数はU ( s)
1= K (1 +)E ( s)
Ti s

(3.3)

ですから

u = K (e +
dz
=e
dt

z
)
Ti

(3.4)

(3.5)

7

オメガシミュレーション

EQUATRAN-Gと直すことができます。

このソーステキストを実行してみると、TREND文によるそして、題意より

積分の途中経過が出力ウィンドウに表示されます。しかし、

e = ys − y(3.6)

数字の羅列では変化の様子がつかみにくいので、グラフに表示したのが図12です。応答の様子がよく分かります。

です。
ソーステキストを図11に示しました。EQUATRAN-Gでは、微分方程式を記述するのに、アポストロフィ（'）を微分記号としてdx
d2x→ x′ ,→ x ′′dt
dt 2のように書き表します。11行目に制御対象の微分方程式(3.2)式が、13～15行目にPI制御器の(3.4)～(3.6)式が入力されています。

図12

例題３の計算結果（グラフ）６．最後にここに示した例題は、説明のためもありますが、式の数が100に満たない簡単な例ばかりでした。実務ではもっと多くの方程式を扱うことになると思いますが、
EQUATRAN-Gは方程式の数が1000を越えるような規模の大きな問題でも解くことができます。
例題３のように、EQUATRAN-Gでダイナミックシミュ
レーションが手軽にできますが、最近はこの技術が応用さ
図11

例題３のソーステキスト

17～18行目は操作量uと弁開度uvの関係を与える条件付きの式で、弁開度の制限
u｜u｜≦ 0.5 のとき｛ 0.5｜u｜＞0.5 のとき

uv＝を表しており、WHEN以降の条件によって場合分けされて計算されます。
23行目のINTEGRAL文が積分の指定をする文で、ｔについて積分区間 0～50で、きざみ 0.05で積分することを指定しています。また、21行目は初期条件の指定で、24、25行
目は出力の指定です。

れてプロセス制御コンピュータなどの実時間システムに組み込まれたり、トレーニングシミュレータのモデル記述用として使われたりなど、その活躍の場をますます広げています。

＜参考文献＞
1)A.L.Myers、W.D.Seider著、大竹伝雄 訳：化学工学の基礎、P.198、培風館、(1982)
2)化学工学協会編：化学工学プログラミング演習、P.154、
培風館、(1976)※EQUATRAN は三井化学（株）の登録商標です。Windows は米国 Microsoft
Corporation の登録商標です。その他の製品名などは一般に各社の商標または登録商標です。

8