fork download 
  1. import numpy as np
  2. from scipy.integrate import solve_ivp
  3. import matplotlib.pyplot as plt
  4.  
  5. # Data
  6.  
  7. # Constants
  8. R = 8.314  # J/mol·K
  9. T0 = 571  # K, reactor inlet temp
  10. T_c0 = 320  # K, coolant inlet temp
  11.  
  12. # Kinetic parameters at 571 K
  13. # k1 = 10.0 /60 # m^3/mol s
  14. # k2 = 5.0  /60 # m^3/mol^0.2 s
  15. # k3 = 2.0  /60 # 1/s
  16.  
  17. # Define Arrhenius function to convert to k(T) in consistent SI units (per s)
  18. def k_arrhenius(T, k0, Ea):
  19.     R = 8.314  # J/mol·K
  20.     return (k0 / 60) * np.exp(-Ea / (R * T))  # convert k0 from per min to per s
  21.  
  22. # Parameters
  23. k0_1, Ea1 = 1.1516e+08, 77.22e3  # m^3/mol·min, J/mol
  24. k0_2, Ea2 = 1.2133e+09, 91.59e3  # mol^(-0.2)·m^(0.6)/min, J/mol
  25. k0_3, Ea3 = 1.2106e+06, 63.32e3  # 1/min, J/mol
  26.  
  27. # Heat of reactions (J/mol)
  28. dH1 = -50000
  29. dH2 = -30000
  30. dHD = 25000
  31.  
  32. # Thermophysical properties
  33. Cp = 150      # J/mol·K (reacting fluid)
  34. rho_c = 1000  # kg/m³
  35. Cp_c = 4180   # J/kg·K (coolant)
  36. U = 100       # W/m²·K
  37.  
  38. # Reactor parameters
  39. v0_total = 100/60    # m^3/s
  40. ca0 = 2              # mol/m^3
  41. n_tubes = 25
  42. d_tube = 20 / 100.0  # m
  43. l_tube = 10          # m
  44. v0 = v0_total / n_tubes  # per tube
  45.  
  46. # Tube geometry
  47. A_tube = np.pi * (d_tube / 2)**2  # m²
  48. V_tube = A_tube * l_tube  # m³
  49. A_surf = np.pi * d_tube * l_tube  # heat transfer area per tube (m²)
  50.  
  51. # Coolant flow rate (kg/s)
  52. m_dot_c = 3.2 /60
  53. m_dot_c = m_dot_c / n_tubes # per tube
  54.  
  55. # Initial conditions: [C_A, C_D, C_U1, C_U2, T, T_coolant]
  56. y0 = [ca0, 0, 0, 0, T0, T_c0]
  57.  
  58. # System of ODEs
  59. def odes(V, y):
  60.  
  61.     C_A, C_D, C_U1, C_U2, T, T_c = y
  62.  
  63.     k1 = k_arrhenius(T, k0_1, Ea1)
  64.     k2 = k_arrhenius(T, k0_2, Ea2)
  65.     k3 = k_arrhenius(T, k0_3, Ea3)
  66.  
  67.     r1 = k1 * C_A**2
  68.     r2 = k2 * C_A**1.2
  69.     r3 = k3 * C_D
  70.  
  71.     # TODO: Write correct expressions for 
  72.     # rates of individual species
  73.     r_A = -r1 -r2
  74.     r_D = 0.5*r1 -r3
  75.     r_U1 = r2
  76.     r_U2 = r3
  77.  
  78.     # TODO: Write correct expressions for 
  79.     # Mole balance
  80.     dC_A_dV = r_A/v0
  81.     dC_D_dV = r_D/v0
  82.     dC_U1_dV = r_U1/v0
  83.     dC_U2_dV = r_U2/v0
  84.  
  85.     # TODO: Write correct expressions for 
  86.     # Energy balance
  87.     Q_gen = (-r1*dH1) + (-r2*dH2) + (r3*dHD)  # J/s
  88.     Q_removal = U * A_surf * (T-T_c)  # W
  89.  
  90.     C_total = C_A + C_D + C_U1 + C_U2 + 1e-6  # avoid zero division
  91.     F_total = C_total * v0  # mol/s
  92.  
  93.     # TODO: Write correct expressions for 
  94.     # Energy balance (reactor side) in terms of Q_gen, Q_removal, and F_total
  95.     dT_dV = (Q_gen - Q_removal) / (m_dot_c * Cp_c) # K/m³
  96.  
  97.     # TODO: Write correct expressions for 
  98.     # Energy balance (coolant side) in terms of Q_removal, m_dot_c, Cp_c
  99.     dT_c_dV = Q_removal / (m_dot_c * Cp_c)  # K/m³
  100.  
  101.     return [dC_A_dV, dC_D_dV, dC_U1_dV, dC_U2_dV,dT_dV,dT_c_dV]
  102.  
  103. # Solve ODEs
  104. V_span = (0, V_tube)
  105. V_eval = np.linspace(V_span[0], V_span[1], 100)
  106.  
  107. solution = solve_ivp(odes, V_span, y0, t_eval=V_eval, method='BDF')
  108. C_A, C_D, C_U1, C_U2, T, T_c = solution.y
  109.  
Success #stdin #stdout #stderr 4.09s 95476KB
comments (?)
stdin 
Standard input is empty
stdout 
Standard output is empty
stderr 
Fontconfig error: No writable cache directories
https://ideone.com/HZC3pt
language:
Python 3 (python  3.12)
created:
2 days ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!