Forums
Blogs
Articles
File Library
FAQ
Search
FB
0
Hi Guys
I am solving these 3 equations to get a temperature and concentration profile in a PFR.
Mass balance: u_(s ) ∂C/∂z=D_(er ) ((∂^2 C)/(∂r^2 )+1/r (∂C_i)/∂r)+D_(er ) ((∂^2 C)/(∂z^2 ))-ρ_c r_BT Energy balance: u_(s ) ρ_f C_(p,f) ∂T/∂z=λ_(er ) ((∂^2 T)/(∂r^2 )+1/r ∂T/∂r)+λ_(er ) ((∂^2 T)/(∂z^2 ))-ΔHρ_c r_BT Reaction rate r_BT=k_0 e^(-E/RT) C
The mass balance includes radius diffusion, axis diffusion and convection and source term and same for energy balance.
I discretise the equation and wrote in Python and it works well. The discretisation form as shown below:
T_(i,j)=1/4 (T_(i+1,j)+T_(i-1,j)+T_(i,j+1)+T_(i,j-1) )+Δr/8r[i] (T_(i+1,j)-T_(i-1,j) )-AΔr/(8λ_(er ) ) (T_(i,j+1)-T_(i,j-1) )-〖Δr〗^2/(4λ_(er ) ) B C_(i,j)=〖1/4(C〗_(i+1,j)+C_(i-1,j)+C_(i,j+1)+C_(i,j-1))+Δr/8r[i] (C_(i+1,j)-C_(i-1,j) )-(u_(s ) Δr)/(8D_(er ) ) (C_(i,j+1)-C_(i,j-1) )-〖〖Δr〗^2/(4D_(er ) ) ρ〗_c 〖*r〗_BT where A = u_(s ) ρ_f C_(p,f) , B= 6×ΔHρ_c r_BT
However, I do not want to include axial diffusion term, so the equatios become like:
Energy Balance: u_(s ) ρ_f C_(p,f) ∂T/∂z=λ_(er ) ((∂^2 T)/(∂r^2 )+1/r ∂T/∂r)-ΔHρ_c r_BT Mass Balance: u_(s ) ∂C/∂z=D_(er ) ((∂^2 C)/(∂r^2 )+1/r (∂C_i)/∂r)-ρ_c r_BT
with discretisation form:
T_(i,j)=1/2 (T_(i+1,j)+T_(i-1,j) )+Δr/4r[i] (T_(i+1,j)-T_(i-1,j) )-AΔr/(4λ_(er ) ) (T_(i,j+1)-T_(i,j-1) )-〖Δr〗^2/(2λ_(er ) ) B C_(i,j)=〖1/2(C〗_(i+1,j)+C_(i-1,j))+Δr/4r[i] (C_(i+1,j)-C_(i-1,j) )-(u_(s ) Δr)/(4D_(er ) ) (C_(i,j+1)-C_(i,j-1) )-〖〖Δr〗^2/(2D_(er ) ) ρ〗_c 〖*r〗_BT
When I delete the axial diffusion term, the code cannot be converged, would you please give some advice? I have been working on it for a week but had no progress..
Code is shown below, it would be much appreciate if you could give some advice..
import numpy as np
import matplotlib.ticker as ticker
import matplotlib.pyplot as plt
from matplotlib.ticker import FuncFormatter
import pandas as pd
from numba import jit, prange
# Set total range
r_range = 0.01 # m
z_range = 0.1 # m
# Set grid size
grid_size = 0.00005 # m
# Calculate the number of points
num_points_r = int(r_range / grid_size) + 1
num_points_z = int(z_range / grid_size) + 1
# Set grid
r_start, r_end = 0, r_range
z_start, z_end = 0, z_range
r = np.linspace(r_start, r_end, num_points_r)
z = np.linspace(z_start, z_end, num_points_z)
# grid distance
dr = r[1] - r[0]
dz = z[1] - z[0]
# initialization of C, T
C = np.ones((num_points_r, num_points_z)) * 3000
T = np.ones((num_points_r, num_points_z)) * 450
Q_reaction = np.zeros((num_points_r, num_points_z))
# Set parameter
u_s = 0.00049 # superficial velocity m/s
Der = 0.000083 # effective radial diffusion coefficient m2/s
rho_catalyst = 4000 # kg/m3
rho_fluid = 704 # kg/m3
cp_fluid = 3037 # J/ kg K
lambda_er = 1.923 # effective radial thermal conductivity, W/ m K
# lambda_er = 6.523 # effective radial thermal conductivity, W/ m K
delta_H = 65000 # J/ mol_H2
Ts = 570 # wall temperature
A = u_s * rho_fluid * cp_fluid
# Set reaction parameter
k = 10816 # m3/kg s
E = 117000 # J/mol
R_gas = 8.314 # J/mol K
# Set initial parameter
C0 = 3622 # mol/m3 rho/MW kg/m3 / kg/mol
T0 = 570 # K
Q0 = k * np.exp(-E / (R_gas * T0)) * C0
Q_reaction[:, :] = Q0 # initial reaction rate, mol_BT/kg_cat s
# Boundary conditions
T[:, 0] = T0 # Inlet
T[:, -1] = T[:, -2] # Outlet
T[0, :] = T[1, :] # lower boundary
T[-1, :] = T0 # constant wall temp, up boundary
C[:, 0] = C0 # Inlet
C[:, -1] = C[:, -2] # Outlet
C[0, :] = C[1, :] # lower boundary
C[-1, :] = C[-2, :] # upper boundary
Q_reaction[:, 0] = Q0
# print("Q_reaction:", Q_reaction)
# iterative to convergence
tolerance = 1e-5
not_converged = True
iteration = 0
max_iterations = 50000000
close_to_convergence = False
# Parallelize code using jit compiler
@jit(nopython=True, parallel=True)
def compute_iteration(C_current, T_current, Q_reaction_current, dr_val, r_val, A_val, lambda_er_val, delta_H_val,
rho_catalyst_val, Der_val, u_s_val, k_val, E_val, R_gas_val):
C_result = np.copy(C_current)
T_result = np.copy(T_current)
Q_reaction_result = np.copy(Q_reaction_current)
for j in prange(1, num_points_z - 1): # loop at z direction
for i in prange(1, num_points_r - 1): # loop at r direction
T_result[i, j] = (1 / 2) * (T_current[i - 1, j] + T_current[i + 1, j]) + (dr_val / (4 * r_val[i])) * \
(T_current[i + 1, j] - T_current[i - 1, j]) - (
(A_val * dr_val) * (T_current[i, j + 1] - T_current[i, j - 1]) / (
4 * lambda_er_val)) + \
(-delta_H_val * 6 * Q_reaction_current[i, j] * rho_catalyst_val * dr_val ** 2) / (
2 * lambda_er_val)
C_result[i, j] = (1 / 2) * (C_current[i - 1, j] + C_current[i + 1, j]) + (dr_val / (4 * r_val[i])) * (
C_current[i + 1, j] - C_current[i - 1, j]) - \
((u_s_val * dr_val) / (4 * Der_val)) * (C_current[i, j + 1] - C_current[i, j - 1]) - (
rho_catalyst_val * Q_reaction_current[i, j]) * \
dr_val ** 2 / (Der_val * 2)
Q_reaction_result[i, j] = k_val * np.exp(-E_val / (R_gas_val * T_result[i, j])) * C_result[i, j]
return C_result, T_result, Q_reaction_result
while not_converged and iteration < max_iterations:
C_new, T_new, Q_reaction_new = compute_iteration(C, T, Q_reaction, dr, r, A, lambda_er, delta_H, rho_catalyst, Der,
u_s, k, E, R_gas)
# Compute the maximum difference between current and previous iterations
max_diff_T = np.max(np.abs(T_new - T))
max_diff_C = np.max(np.abs(C_new - C))
if (iteration + 1) % 1000 == 0:
print(
f"Iteration {iteration + 1}: Max difference in T = {max_diff_T:.10f}, Max difference in C = {max_diff_C:.10f}")
# If the maximum difference is below our threshold, set the flag to check for convergence in the next iteration
if max_diff_T < tolerance and max_diff_C < tolerance:
close_to_convergence = True
# Check for convergence if close_to_convergence is True
if close_to_convergence:
not_converged = max_diff_T > tolerance or max_diff_C > tolerance
# Update C, T, and Q_reaction for the next iteration
C, T, Q_reaction = C_new, T_new, Q_reaction_new
# Boundary conditions remain the same
T[:, 0] = T0 # Inlet
T[:, -1] = T[:, -2] # Outlet
T[0, :] = T[1, :] # lower boundary
T[-1, :] = T0 # constant wall temp, up boundary
C[:, 0] = C0 # Inlet
C[:, -1] = C[:, -2] # Outlet
C[0, :] = C[1, :] # lower boundary
C[-1, :] = C[-2, :] # upper boundary
Q_reaction[:, 0] = Q0
iteration += 1
fig, ax = plt.subplots(2, 1, figsize=(10, 10)) # 3 rows, 1 column
# Temperature plot
img1 = ax[0].imshow(T, cmap='hot', interpolation='nearest', origin='lower', vmin=np.min(T), vmax=np.max(T),
aspect=0.5 * len(z) / len(r), extent=[z_start, z_end, r_start, r_end])
fig.colorbar(img1, ax=ax[0], format=ticker.FuncFormatter(lambda x, pos: '{:.10f}'.format(x)), label='Temperature (K)')
ax[0].set_title('Heat Distribution')
ax[0].set_xlabel('z')
ax[0].set_ylabel('r')
# Concentration plot
img2 = ax[1].imshow(C, cmap='hot', interpolation='nearest', origin='lower', vmin=np.min(C), vmax=np.max(C),
aspect=0.5 * len(z) / len(r), extent=[z_start, z_end, r_start, r_end])
fig.colorbar(img2, ax=ax[1], format=ticker.FuncFormatter(lambda x, pos: '{:.10f}'.format(x)),
label='Concentration (mol/m3)')
ax[1].set_title('Concentration')
ax[1].set_xlabel('z')
ax[1].set_ylabel('r')
plt.tight_layout()
plt.show()
r_index_half = num_points_r // 2 # Index corresponding to r = r_range/2
plt.figure(figsize=(12, 6))
# Subplot for Temperature
plt.subplot(1, 2, 1) # 1 row, 2 columns, first subplot
plt.plot(z, T[0, :], label=f'r = 0')
plt.plot(z, T[r_index_half, :], label=f'r = r_range/2')
plt.plot(z, T[-1, :], label=f'r = r_range')
plt.title("Temperature along z axis")
plt.xlabel("z")
plt.gca().yaxis.set_major_formatter(FuncFormatter(lambda x, pos: '{:.10f}'.format(x))) # Change to decimal format
plt.ylabel("Temperature")
plt.legend()
# Subplot for Concentration
plt.subplot(1, 2, 2) # 1 row, 2 columns, second subplot
plt.plot(z, C[0, :], label=f'r = 0')
plt.plot(z, C[r_index_half, :], label=f'r = r_range/2')
plt.plot(z, C[-1, :], label=f'r = r_range')
plt.title("Concentration along z axis")
plt.xlabel("z")
plt.gca().yaxis.set_major_formatter(FuncFormatter(lambda x, pos: '{:.10f}'.format(x))) # Change to decimal format
plt.ylabel("Concentration")
plt.legend()
plt.tight_layout()
plt.show()
# Plot 4
plt.plot(z, Q_reaction[r_index_half, :], label='MCH')
plt.title("Reaction rate")
plt.xlabel("z")
plt.gca().yaxis.set_major_formatter(FuncFormatter(lambda x, pos: '{:.10f}'.format(x))) # Change to decimal format
plt.ylabel("Rate")
plt.legend()
plt.show()
