ExcelマクロVBAで回分式反応器の設計計算~反応速度・物質収支・熱収支を組み込んだ連立微分方程式の解析~





ExcelマクロVBAで回分式反応器の設計計算~反応速度・物質収支・熱収支を組み込んだ連立微分方程式の解析~



資料請求番号:SH43 TS56

反応器設計の基礎的な計算を行うプログラムをVBAで作成

今回の資料は↓こちら↓の続きで、↓のページで導いた、回分式反応器の濃度と温度に関する連立微分方程式をVBAで解いて、課題を解決するということを趣旨としています。
[blogcard url=”https://shimaphoto03.com/science/batch/”]

課題① A→Rタイプ

『A→Rで表される1次反応を、下記に示すような回分反応器で液相中で行う。
反応器は外側ジャケットに393Kのオイルを流して保温する。反応は発熱反応である。この時、反応物濃度と槽内温度の経時変化を求めよ。』

この反応系の基礎式は



です。
基礎式はこちらで導出しています。

従って、反応容積、反応器/ジャケット間の総括熱伝達係数、伝熱面積といった、設備の情報
反応液密度、反応液比熱、反応の頻度因子と活性化エネルギーといった、物質の情報
微分方程式を解くための初期条件として、初期濃度とジャケットに流入するオイルの温度といった、操作情報

が必要になります。今回は、以下のように設定しました。また、反応容積、伝熱面積については直径50cm, 高さ1mの反応器を使用することとし、液位は高さの70%にすることとしました。

マクロを使用するExcelシートスタイル

マクロを適用するExcelシートのスタイルは以下の通りです。単位はすべてSI基本単位系に従っています。

プログラム

上記のスタイルで下記のマクロを実行すると上記連立微分方程式の解である濃度と温度の経時変化を取得できます。

※予めグラフの範囲の指定をしておいてください。

Public k0, E, R, U, A, rho, V, Cp, Tw, deltaH, rA As Double

Sub batch()

Dim CA0, Temp0 As Double
Dim rA As Double
Dim init, ed, h, h2, kCA(4), kTemp(4) As Double
Dim i, j As Integer

‘CalculaTemp propertv
init = Cells(1, 3)
ed = Cells(2, 3)
h = Cells(3, 3)
h2 = Cells(4, 3)

‘Constants
V = Cells(5, 3)
U = Cells(6, 3)
A = Cells(7, 3)
rho = Cells(8, 3)
Cp = Cells(9, 3)
k0 = Cells(10, 3)
E = Cells(11, 3)
deltaH = Cells(12, 3)
CA0 = Cells(13, 3)
Tw = Cells(14, 3)
R = Cells(15, 3)

‘Initial Condition
CA0 = Cells(4, 7)
Temp0 = Cells(4, 8)
CA = CA0
Temp = Temp0
t = init

‘Runge-kutta roop

For i = 0 To ((ed – init) / h)

j = 6 + i / h2

Cells(1, 6) = i
Cells(2, 6) = j – 6

If i Mod h2 = 0 Then

Cells(j, 6) = t
Cells(j, 7) = CA
Cells(j, 8) = Temp
Cells(j, 9) = k0 * Exp(-1 * E / (R * Temp)) * CA

End If

kCA(1) = h * F1(t, CA, Temp)
kTemp(1) = h * F2(t, CA, Temp)
kCA(2) = h * F1(t + h / 2, CA + kCA(1) / 2, Temp + kTemp(1) / 2)
kTemp(2) = h * F2(t + h / 2, CA + kCA(1) / 2, Temp + kTemp(1) / 2)
kCA(3) = h * F1(t + h / 2, CA + kCA(2) / 2, Temp + kTemp(2) / 2)
kTemp(3) = h * F2(t + h / 2, CA + kCA(2) / 2, Temp + kTemp(2) / 2)
kCA(4) = h * F1(t + h, CA + kCA(3), Temp + kTemp(3))
kTemp(4) = h * F2(t + h, CA + kCA(3), Temp + kTemp(3))
nCA = CA + (kCA(1) + 2 * kCA(2) + 2 * kCA(3) + kCA(4)) / 6
nTemp = Temp + (kTemp(1) + 2 * kTemp(2) + 2 * kTemp(3) + kTemp(4)) / 6
nt = t + h

t = nt
CA = nCA
Temp = nTemp

Next i

End Sub

Function F1(ByVal t As Double, ByVal CA As Double, ByVal Temp As Double) As Double
rA = k0 * Exp(-1 * E / (R * Temp)) * CA
F1 = -rA
End Function

Function F2(ByVal t As Double, ByVal CA As Double, ByVal Temp As Double) As Double
F2 = -(U * A / (rho * V * Cp)) * (Temp – Tw) + (-deltaH / (rho * Cp)) * rA
End Function

実行後の結果

課題① 応用例

今回の計算により、393K=120℃のオイルを用いて運転すれば約5,000秒(=1時間20分)で50%の、16900秒(=4時間41分)で80%の反応率が得られることが分かった。温度上昇も2℃以下であり、熱暴走の危険もほとんどない。

ところが、反応に時間がかかりすぎていることが問題視され、

「できるだけ早く反応を済ませることはできないか?」

という要望が出された。本マクロを使用して、反応条件の再提案を行え。

回答例

早く反応を済ませるのに一番最初に思いつきうる方法として
「温度を上げて運転すること」であるが、温度を上げることによって

・反応器の耐熱能力を上回った状態で操作してしまう
・反応物質の蒸気圧が反応器の耐圧を超えてしまう

というリスクが考えられたので、「反応器の耐熱温度」「反応器の耐圧性能」を問合せ、「反応物質の蒸気圧」を求めたところ、
473K(=200℃)までであれば、安全に運転できるが、余裕を見て453K(=180℃)まで上げてよいという条件を設定した。

そこで、以下のマクロを実行して393K~453Kまでシミュレーションをおこなった。

Sub batchstep()

For Tw = 393 To 453 Step 10
Cells(14, 3) = Tw
Application.Run “batch”
DoEvents: DoEvents: DoEvents
Next Tw

End Sub

このシミュレーション結果から423K~443Kまでの保温温度で詳細に反応器内温度、濃度を調べることにした。積分範囲を0~20000sから0~3000sまでと改め、刻み幅も10sから1sへ変更した。

この動画より434Kまで保温温度を上げてもよいことが分かった。この温度で運転する場合、340s=(5分40秒)で反応率90%を達成できることがわかる。

課題② A→R→Sタイプ

『A→R→Sで表される1次の逐次反応を、下記に示すような回分反応器で液相中で行う。
反応器は外側ジャケットに393Kのオイルを流して保温する。反応は発熱反応である。この時、反応物濃度と槽内温度の経時変化を求めよ。』

この反応系の基礎式は




です。

基礎式はこちらで導出しています。

前回のA→R型反応と比較して、反応の数が増えたので、それに対応するように頻度因子と活性化エネルギーと反応エンタルピーが、
解析したい物質の個数が増えたので、各物質の初期濃度の情報が増えます。

ただし、今回は反応液の密度、比熱、容積が一定であると仮定を置いているため、これらの変数の数は変化しません。

シートスタイルとプログラム

資料が冗長になってしまうので、たたんでおきます。
”シートスタイル”



”プログラム”


Public k10, k20, E1, E2, R, U, A, rho, V, Cp, Tw, dH1, dH2, r1, r2 As Double

Sub batch2()

Dim CA0, CR0, CS0, Temp0 As Double
Dim CA, CR, CS, Temp As Double
Dim init, ed, h, h2, kCA(4), kCR(4), kCS(4), kTemp(4) As Double
Dim i, j As Integer

‘CalculaTemp propertv
init = Cells(1, 3)
ed = Cells(2, 3)
h = Cells(3, 3)
h2 = Cells(4, 3)

‘Constants
V = Cells(5, 3)
U = Cells(6, 3)
A = Cells(7, 3)
rho = Cells(8, 3)
Cp = Cells(9, 3)
k10 = Cells(10, 3)
k20 = Cells(11, 3)
E1 = Cells(12, 3)
E2 = Cells(13, 3)
dH1 = Cells(14, 3)
dH2 = Cells(15, 3)
CA0 = Cells(16, 3)
CR0 = Cells(17, 3)
CS0 = Cells(18, 3)
Tw = Cells(19, 3)
R = Cells(20, 3)

‘Initial Condition
CA0 = Cells(4, 7)
CR0 = Cells(4, 8)
CS0 = Cells(4, 9)
Temp0 = Cells(4, 10)

CA = CA0
CR = CR0
CS = CS0
Temp = Temp0
t = init

‘Runge-kutta roop

For i = 0 To ((ed – init) / h)

j = 6 + i / h2

Cells(1, 6) = i
Cells(2, 6) = j – 6

If i Mod h2 = 0 Then

Cells(j, 6) = t
Cells(j, 7) = CA
Cells(j, 8) = CR
Cells(j, 9) = CS
Cells(j, 10) = Temp
Cells(j, 11) = kCS(1)
Cells(j, 12) = CR

End If

kCA(1) = h * F1(t, CA, CR, CS, Temp)
kCR(1) = h * F2(t, CA, CR, CS, Temp)
kCS(1) = h * F3(t, CA, CR, CS, Temp)
kTemp(1) = h * F4(t, CA, CR, CS, Temp)

kCA(2) = h * F1(t + h / 2, CA + kCA(1) / 2, CR + kCR(1) / 2, CS + kCS(1) / 2, Temp + kTemp(1) / 2)
kCR(2) = h * F2(t + h / 2, CA + kCA(1) / 2, CR + kCR(1) / 2, CS + kCS(1) / 2, Temp + kTemp(1) / 2)
kCS(2) = h * F3(t + h / 2, CA + kCA(1) / 2, CR + kCR(1) / 2, CS + kCS(1) / 2, Temp + kTemp(1) / 2)
kTemp(2) = h * F4(t + h / 2, CA + kCA(1) / 2, CR + kCR(1) / 2, CS + kCS(1) / 2, Temp + kTemp(1) / 2)

kCA(3) = h * F1(t + h / 2, CA + kCA(2) / 2, CR + kCR(2) / 2, CS + kCS(2) / 2, Temp + kTemp(2) / 2)
kCR(3) = h * F2(t + h / 2, CA + kCA(2) / 2, CR + kCR(2) / 2, CS + kCS(2) / 2, Temp + kTemp(2) / 2)
kCS(3) = h * F3(t + h / 2, CA + kCA(2) / 2, CR + kCR(2) / 2, CS + kCS(2) / 2, Temp + kTemp(2) / 2)
kTemp(3) = h * F4(t + h / 2, CA + kCA(2) / 2, CR + kCR(2) / 2, CS + kCS(2) / 2, Temp + kTemp(2) / 2)

kCA(4) = h * F1(t + h, CA + kCA(3), CR + kCR(3), CS + kCS(3), Temp + kTemp(3))
kCR(4) = h * F2(t + h, CA + kCA(3), CR + kCR(3), CS + kCS(3), Temp + kTemp(3))
kCS(4) = h * F3(t + h, CA + kCA(3), CR + kCR(3), CS + kCS(3), Temp + kTemp(3))
kTemp(4) = h * F4(t + h, CA + kCA(3), CR + kCR(3), CS + kCS(3), Temp + kTemp(3))

nCA = CA + (kCA(1) + 2 * kCA(2) + 2 * kCA(3) + kCA(4)) / 6
nCR = CR + (kCR(1) + 2 * kCR(2) + 2 * kCR(3) + kCR(4)) / 6
nCS = CS + (kCS(1) + 2 * kCS(2) + 2 * kCS(3) + kCS(4)) / 6
nTemp = Temp + (kTemp(1) + 2 * kTemp(2) + 2 * kTemp(3) + kTemp(4)) / 6
nt = t + h

t = nt
CA = nCA
CR = nCR
CS = nCS
Temp = nTemp

Next i

End Sub

Function F1(ByVal t As Double, ByVal CA As Double, ByVal CR As Double, ByVal CS As Double, ByVal Temp As Double) As Double
r1 = k10 * Exp(-1 * E1 / (R * Temp)) * CA
r2 = k20 * Exp(-1 * E2 / (R * Temp)) * CR
F1 = -r1
End Function

Function F2(ByVal t As Double, ByVal CA As Double, ByVal CR As Double, ByVal CS As Double, ByVal Temp As Double) As Double
F2 = r1 – r2
End Function

Function F3(ByVal t As Double, ByVal CA As Double, ByVal CR As Double, ByVal CS As Double, ByVal Temp As Double) As Double
F3 = r2
End Function

Function F4(ByVal t As Double, ByVal CA As Double, ByVal CR As Double, ByVal CS As Double, ByVal Temp As Double) As Double
F4 = -(U * A / (rho * V * Cp)) * (Temp – Tw) + (-dH1 / (rho * Cp)) * r1 + (-dH2 / (rho * Cp)) * r2
End Function

実行の結果

まとめ

今回は、回分反応器設計の基礎的な計算を行うプログラムをVBAで作成し、回分反応器に関する「反応速度式」「物質収支式」「熱収支式」をVBAで解いた結果を示しました。

shimakei8364