(k :反応速度定数, p :A,B成分分圧[kPa], R :気体定数 8.314J/(mol･K), T :温度[K])
（久保田・関沢「反応工学概論」P.149,例題12.2）
量論式A+B→Rなる気相触媒反応の速度 r
は次のように与えられる。
(k :反応速度定数, p
:A,B成分分圧[kPa], R
:気体定数 8.314J/(mol･K), T :温度[K])
A成分の転化率x と各成分分圧の関係は，
(P :全圧[kPa], zA0,zB0:A,B成分入り口モル分率,δA:量論係数(=(1-1-1)/1)
)
全圧 P =101.3kPa の条件でAが 5%, Bが 20%その他は不活性成分の気体原料を,流量F =(35000+100*(学籍番号下2桁))mol/h ,温度673Kで反応装置に供給する。原料Aの転化率 x [-]および供給気体の温度 T [K]の流れ方向の変化を記述する微分方程式は次式となる。
(ρB=1100kg/m3:触媒密度, Q=-142.4kJ/mol:反応熱, Cp=31.4J/(mol･K)：供給気体熱容量, U =2.5x105J/(m2･h･K):総括伝熱係数, A
:単位長さ当り伝熱面積[m], S
:流れ方向断面積[m2], Tw:冷却媒体温度）
T の初期値673Kから出発して (UA/S) = 0 , 5x106, 5x107 の条件で微分方程式を積分し，横軸 V [m3](0～2)対縦軸 T [K]および x のグラフを描く。
① 最高温度を750K以下におさえる条件を探し，そのときx >99%の条件で触媒層容積 V
[m3]を決める。
②(UA/S)
の値と S/A
= n (1/4)πd 2
/ nπd =
(1/4)d
の関係から触媒管径d [m]を求める。
③長さ l [m]を仮定してV =n(1/4)πd 2l より触媒管本数n が求まる。反応管を2d ピッチ三角配列とすると，容器径D [m]の反応容器中で一本の反応管が占める面積はπd 2と近似できるので，nπd 2=(1/4)πD 2 よりD が決る。なるべくl ≒D となるようにD, l
を決める。
変数の割り当て：転化率 x → Y(0) , 温度 T
→Y(1)
1000 '====ﾚﾝﾘﾂ ｼﾞｮｳﾋﾞﾌﾞﾝ
ﾎｳﾃｲｼｷ ========
1010 ' RUNGE-KUTTA METHOD FOR SYSTEMS OF DIFFERENTIAL EQUATIONS
1020 READ NM
1030 NM = NM - 1
1040 DIM Y(NM), Y0(NM), YW(NM), YS(NM), YQ(NM), YD(NM)
1050 READ YA, YB
1070 FOR NK = 0 TO NM
中間略
2000 '---------ﾋﾞﾌﾞﾝ ﾎｳﾃｲｼｷ
2010 '----ﾎｳﾃｲｼｷ ﾉ ｶｽﾞ
2020 DATA 2
2030 '----ｾｷﾌﾞﾝ ｸｶﾝ [a,b]
2040 DATA 0,2.0
2050 '----ｼｮｷﾁ x=a ﾆｵｹﾙ Y(0),Y(1),Y(2)...ﾉ
ｱﾀｲ
2060 DATA 0.0 , 673
2070 '----xﾉ ｷｻﾞﾐﾊﾊﾞh, ｺﾞｻ ﾉ ｷｮﾖｳ
ｹﾞﾝｶｲ E
2080 DATA 0.1 , 1
2085 '----EQUATIONS( Y,X,N ﾃﾞﾊｼﾞﾏﾙ ﾍﾝｽｳ ﾊ
ﾂｶﾜﾅｲ)
2090 F0 = 40000 : UAS= 0
各自修正
2095 ZA0=.05 ZB0=.2 : P=101.3 : DELTAA=(1-1-1)/1
2096 CP=31.4 : Q=-142400! : TW=673 : ROB=1100
2097 PA = ZA0*(1-Y(0))*P /(1+DELTAA*ZA0*Y(0))
2098 PB = (ZB0-ZA0*Y(0))*P /(1+DELTAA*ZA0*Y(0))
2099 K = 24400!*EXP(-75370!/8.314/Y(1)) : R=-K*PA*PB
2100 ON NK + 1 GOTO 2110, 2130
2110 ' [ DYDX=Y(0)' ]
2120 DYDX =- R*ROB /(ZA0*F0)
2125 RETURN
2130 ' [ DYDX=Y(1)' ]
2140 DYDX = ( Q*R*ROB-UAS*(Y(1)-TW)) / ( CP*F0)
2145 RETURN