資料請求番号:SH43 TS56
今回の資料は↓こちら↓の続きで、↓のページで導いた、回分式反応器の濃度と温度に関する連立微分方程式をVBAで解いて、課題を解決するということを趣旨としています。
[blogcard url=”https://shimaphoto03.com/science/batch/”]
『A→Rで表される1次反応を、下記に示すような回分反応器で液相中で行う。
反応器は外側ジャケットに393Kのオイルを流して保温する。反応は発熱反応である。この時、反応物濃度と槽内温度の経時変化を求めよ。』
この反応系の基礎式は
です。
基礎式はこちらで導出しています。
従って、反応容積、反応器/ジャケット間の総括熱伝達係数、伝熱面積といった、設備の情報、
反応液密度、反応液比熱、反応の頻度因子と活性化エネルギーといった、物質の情報、
微分方程式を解くための初期条件として、初期濃度とジャケットに流入するオイルの温度といった、操作情報
が必要になります。今回は、以下のように設定しました。また、反応容積、伝熱面積については直径50cm, 高さ1mの反応器を使用することとし、液位は高さの70%にすることとしました。
マクロを適用する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で表される1次の逐次反応を、下記に示すような回分反応器で液相中で行う。
反応器は外側ジャケットに393Kのオイルを流して保温する。反応は発熱反応である。この時、反応物濃度と槽内温度の経時変化を求めよ。』
この反応系の基礎式は
です。
基礎式はこちらで導出しています。
前回のA→R型反応と比較して、反応の数が増えたので、それに対応するように頻度因子と活性化エネルギーと反応エンタルピーが、
解析したい物質の個数が増えたので、各物質の初期濃度の情報が増えます。
ただし、今回は反応液の密度、比熱、容積が一定であると仮定を置いているため、これらの変数の数は変化しません。
資料が冗長になってしまうので、たたんでおきます。 ”プログラム” Sub batch2() Dim CA0, CR0, CS0, Temp0 As Double ‘CalculaTemp propertv ‘Constants ‘Initial Condition CA = CA0 ‘Runge-kutta roop For i = 0 To ((ed – init) / h) j = 6 + i / h2 Cells(1, 6) = i If i Mod h2 = 0 Then Cells(j, 6) = t End If kCA(1) = h * F1(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) kCA(3) = h * F1(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)) nCA = CA + (kCA(1) + 2 * kCA(2) + 2 * kCA(3) + kCA(4)) / 6 t = nt 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 Function F2(ByVal t As Double, ByVal CA As Double, ByVal CR As Double, ByVal CS As Double, ByVal Temp As Double) As Double Function F3(ByVal t As Double, ByVal CA As Double, ByVal CR As Double, ByVal CS As Double, ByVal Temp As Double) As Double Function F4(ByVal t As Double, ByVal CA As Double, ByVal CR As Double, ByVal CS As Double, ByVal Temp As Double) As Double 実行の結果 今回は、回分反応器設計の基礎的な計算を行うプログラムをVBAで作成し、回分反応器に関する「反応速度式」「物質収支式」「熱収支式」をVBAで解いた結果を示しました。
”シートスタイル”
Public k10, k20, E1, E2, R, U, A, rho, V, Cp, Tw, dH1, dH2, r1, r2 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
init = Cells(1, 3)
ed = Cells(2, 3)
h = Cells(3, 3)
h2 = Cells(4, 3)
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)
CA0 = Cells(4, 7)
CR0 = Cells(4, 8)
CS0 = Cells(4, 9)
Temp0 = Cells(4, 10)
CR = CR0
CS = CS0
Temp = Temp0
t = init
Cells(2, 6) = j – 6
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
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)
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)
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)
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))
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
CA = nCA
CR = nCR
CS = nCS
Temp = nTemp
r1 = k10 * Exp(-1 * E1 / (R * Temp)) * CA
r2 = k20 * Exp(-1 * E2 / (R * Temp)) * CR
F1 = -r1
End Function
F2 = r1 – r2
End Function
F3 = r2
End Function
F4 = -(U * A / (rho * V * Cp)) * (Temp – Tw) + (-dH1 / (rho * Cp)) * r1 + (-dH2 / (rho * Cp)) * r2
End Functionまとめ