🐙

Julia : ModelingToolkit での 連続撹拌槽型反応器を解析

に公開
2025/12/31
Julia
化学工学
modelingtoolkit
反応工学
tech

Julia の ModelingToolkit を使って、連続撹拌槽型反応器の反応パフォーマンスを求める問題を解いていました。
問題自体は、簡単な代数方程式になりますが、実際の反応装置、条件に合わせてパラメメータや変数を入れ替えて計算するのに便利かなとおもい整理しました。

定常状態の完全混合連続撹拌槽型反応器の支配方程式

体積（密度）変化のない液相での反応系として、完全混合連続撹拌槽型の設計方程式について、成分Aに着目した場合には、

input = output + disapperance of reaction + accumulation

input : v C_{A,0}

output : v C_A

disapperance of reaction : ( -r_A )V

accumulation : V \frac{d C_A}{dt}

今回は、定常状態の計算を行うので、accumulation（時間依存）の項は、ゼロになります。

v C_{A,0} = v C_A + (-r_A) V + 0

となります。
v : 体積流量 [L/min]
V : 反応器体積 [L]
C_A : 反応器内の成分Aの濃度 [mol/L]
C_{A,0} : 成分Aの供給濃度
r_A : 反応速度式 [mol/min-L]

例題

\text{A + 2B} \rightleftharpoons \text{R}

成分Aに対する反応速度式を以下で表せるとする。

-r_A = k_1 C_A C_B^2 - k_2 C_R

k_1 : 正反応の反応速度定数 [L^2/min-mol^2]
k_2 : 逆反応の反応速度定数 [1/min]

とする。

操作条件を以下とする。
反応器体積 V= 6 [m^3]
各原料A,Bを v = 1.0 [L/min]
初期濃度 C_{A,0} = 2.8 [mol/L]、C_{B,0} = 1.6 [mol/L]、C_{R,0} = 0 [mol/L]
反応速度定数 k_1 = 12.5 [L^2/min-mol^2]、k_2 = 1.5 [1/min]
とする。

この条件での出口の濃度を求める。関係式を整理すると

-r_A = -\frac{1}{2} r_B
-r_A = r_R
-r_A = k_1 C_A C_B^2 - k_2 C_R
v C_{A,0} = v C_A + ( -r_A ) V
v C_{B,0} = v C_B + ( -r_B ) V
v C_{R,0} = v C_R + ( -r_R ) V
\tau = \frac{V}{v}
X_A = 1 - \frac{C_A}{C_{A,0}}
X_B = 1 - \frac{C_A}{C_{B,0}}

となる。

ここで、変数(variables)は、
r_A, r_B, r_R, C_A, C_B, C_R, \tau, X_A, X_B

与条件（parameters)
k_1, k_2, V, v, C_{A,0}, C_{B,0}, C_{R,0}
である。

ModelingToolkit での解放

供給体積流量をv0 としている。A,Bを1 L/min で供給するので、2 L/min となる。また初期濃度は、体積が倍になるので、A,Bの組成はその半分でスタートするとこになる。

v0 = 2.0 [L/min]
C_{A,0} = 1.4 [mol/L]
C_{B,0} = 0.8 [mol/L]

となる。

using ModelingToolkit, NonlinearSolve

# Define a nonlinear system
@variables CA CB CR rA rB rR τ XA XB
@parameters k1 k2 V CA0 CB0 CR0 v0
eqs = [rA ~ 1/2 * rB
       rA ~ - k1 * CA * CB^2 + k2 * CR
       rA ~ -rR
       v0 * CA0 ~ v0 * CA + (-rA) * V
       v0 * CB0 ~ v0 * CB + (-rB) * V
       v0 * CR0 ~ v0 * CR + (-rR) * V
       τ ~ V / v0 #(CA0 - CA) / (-rA)
       XA ~ 1 - CA / CA0
       XB ~ 1 - CB / CB0]
@mtkcompile ns = System(eqs)

guesses = [CA => 0.1, CB => 0.1, CR => 1.0, rA => 0.1, rB => 0.2, rR => 0.1, τ => 1.0, XA => 1.0]
ps = [k1 => 12.5, k2 => 1.5, V => 6, CA0 => 1.4, CB0 => 0.8, CR0 => 0, v0 => 2.0]   

prob = NonlinearProblem(ns, vcat(guesses, ps))
sol = solve(prob, NewtonRaphson())
答え
retcode: Success
u: 3-element Vector{Float64}:
 1.0999999999999999
 0.20000000000000004
 0.1

独立変数を確認すると

unknowns(ns)
答え
3-element Vector{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}}:
 CA
 CB
 rR

以上の３つが独立変数で他の変数は従属変数としてもとまる。

println("CA = ", sol[CA])
println("CB = ", sol[CB])
println("CR = ", sol[CR])
println("rA = ", sol[rA])
println("rB = ", sol[rB])
println("rR = ", sol[rR])
println("τ = ", sol[τ])
println("XA = ", sol[XA])
println("XB = ", sol[XB])
答え
CA = 1.0999999999999999
CB = 0.20000000000000004
CR = 0.30000000000000004
rA = -0.10000000000000009
rB = -0.20000000000000018
rR = 0.1
τ = 3.0
XA = 0.2142857142857143
XB = 0.75

例題２

先の問題では、流量を与えて、Bの転化率を求める計算を実施した。
今度は、Bの転化率を0.75として流量を求める問題としてみる。

XB を parameter として、0.75 を与えて、 v0 を variables として適当な初期値を与えて計算する。

プログラムは以下になる。

using ModelingToolkit, NonlinearSolve

# Define a nonlinear system
@variables CA CB CR rA rB rR τ XA v0
@parameters k1 k2 V CA0 CB0 CR0 XB
eqs = [rA ~ 1/2 * rB
       rA ~ - k1 * CA * CB^2 + k2 * CR
       rA ~ -rR
       v0 * CA0 ~ v0 * CA + (-rA) * V
       v0 * CB0 ~ v0 * CB + (-rB) * V
       v0 * CR0 ~ v0 * CR + (-rR) * V
       τ ~ V / v0 #(CA0 - CA) / (-rA)
       XA ~ 1 - CA / CA0
       XB ~ 1 - CB / CB0
       ]
@mtkcompile ns = System(eqs)

guesses = [CA => 0.1, CB => 0.1, CR => 1.0, rA => 0.1, rB => 0.2, rR => 0.1, τ => 1.0, XA => 1.0, v0 => 2.0]
ps = [k1 => 12.5, k2 => 1.5, V => 6, CA0 => 1.4, CB0 => 0.8, CR0 => 0, XB => 0.75]   

prob = NonlinearProblem(ns, vcat(guesses, ps))
sol = solve(prob, NewtonRaphson())
答え
retcode: Success
u: 2-element Vector{Float64}:
 1.100000000000001
 2.0000000000000018

独立変数を確認する。

unknowns(ns)
答え
2-element Vector{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}}:
 CA
 v0

他の変数を求める。XBはparameterであるので、コメントとして、v0 を追加している。

println("CA = ", sol[CA])
println("CB = ", sol[CB])
println("CR = ", sol[CR])
println("rA = ", sol[rA])
println("rB = ", sol[rB])
println("rR = ", sol[rR])
println("τ = ", sol[τ])
println("XA = ", sol[XA])
# println("XB = ", sol[XB])
println("v0 = ", sol[v0])
答え
CA = 1.100000000000001
CB = 0.2
CR = 0.30000000000000054
rA = -0.0999999999999998
rB = -0.20000000000000018
rR = 0.10000000000000027
τ = 2.9999999999999973
XA = 0.21428571428571352
v0 = 2.0000000000000018

まとめ

化学工学でよく使われる反応器のモデルの一つである連続撹拌槽型反応器を例にMoleingToolkitを利用した。

化学工学の背景、理論の詳細は割愛したので、追って整理をする。

実行環境

versioninfo()

Julia Version 1.12.3
Commit 966d0af0fdf (2025-12-15 11:20 UTC)
Build Info:
  Official https://julialang.org release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 8 × AMD Ryzen 5 3500C with Radeon Vega Mobile Gfx
  WORD_SIZE: 64
  LLVM: libLLVM-18.1.7 (ORCJIT, znver1)
  GC: Built with stock GC
Threads: 1 default, 1 interactive, 1 GC (on 8 virtual cores)
Environment:
  LD_ARGV0_REL = ../bin/garcon

Discussion