求指导 可以加钱！
固定床催化反应，二维稳态传热传质问题，有源项，解线性偏微分方程
请问我这个离散做得对吗
传热方程 u_s * rho_fluid * Cp_fluid * ∂T/∂z = λ_effective * (∂^2T/∂r^2 + 1/r * ∂T/∂r) + (-ΔH)* r * ρ_catalyst
传质方程 u_s * ∂C/∂z = D_effective * (∂^2C/∂r^2 + 1/r * ∂C/∂r) + ρ_catalyst * r
其中r是反应速率方程 一阶反应 r = k0*e^(-E/RT) * C
r为径向 z为轴向 一个固定床催化反应
我用中心差分一直收敛不了 不知道哪里出错了
import numpy as np
import matplotlib.pyplot as plt
# Set grid
r_start, r_end, num_points_r = 0, 100, 100
z_start, z_end, num_points_z = 0, 100, 100
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.zeros((num_points_r, num_points_z))
T = np.zeros((num_points_r, num_points_z))
# Initialize lists to store changes
residuals_C = []
residuals_T = []
# Set parameter
u_s = 0.0049 # superficial velocity m/s
Der = 0.00000895 # effective radial diffusion coefficient m2/s
rho_catalyst = 4000 # kg/m3
rho_fluid = 704 # kg/m3
cp_fluid = 3037 # J/ kg K
lambda_er = 3.98 # effective radial thermal conductivity, W/ m K
delta_H = -65000 # J/ mol_H2
Ut = 10 # heat transfer coefficient W/ m2 K
Ts = 570 # wall temperature
# Set reaction parameter
k = 0.0108 # m3/kg s
E = 117000 # J/mol
R_gas = 8.314 # J/mol K
# Set initial parameter
C0 = 3694 # mol/m3 rho/MW
T0 = 570 # K
C[:, 0] = C0
T[:, 1] = T0
# iterative to convergence
tolerance = 1e-9
not_converged = True
iteration = 0
max_iterations = 1000
reaction_i = k * np.exp(-E / (R_gas * Ts))
print("1st ri", reaction_i)
while not_converged and iteration < max_iterations:
# Initialize C_new and T_new with the current state of C and T
C_new = np.copy(C)
T_new = np.copy(T)
for j in range(1, num_points_z - 1): # loop at z direction
for i in range(1, num_points_r - 1): # loop at r direction
# Calculate reaction rate
reaction_i = k * np.exp(-E / (R_gas * T[i, j])) # m3/kg s
print("T", reaction_i)
# Central Finite Difference Approximations, r, 2 order derivative
d2C_dr2 = (C[i + 1, j] - 2 * C[i, j] + C[i - 1, j]) / dr ** 2
d2T_dr2 = (T[i + 1, j] - 2 * T[i, j] + T[i - 1, j]) / dr ** 2
# Central Finite Difference Approximations, r, 1 order derivative
dC_dr = (C[i + 1, j] - C[i - 1, j]) / (2 * dr)
dT_dr = (T[i + 1, j] - T[i - 1, j]) / (2 * dr)
# Central Finite Difference Approximations, z, 1 order derivative
dC_dz = (C[i, j + 1] - C[i, j - 1]) / (2 * dz)
dT_dz = (T[i, j + 1] - T[i, j - 1]) / (2 * dz)
# update C_new[i, j]
C_new[i, j] = (u_s * dC_dz - Der * (d2C_dr2 + (1 / r[i]) * dC_dr)) / (
rho_catalyst * reaction_i)
# update T_new[i, j]
T_new[i, j] = (1 / 2) * (
((u_s * rho_fluid * cp_fluid * (dr ** 2)) / lambda_er) * dT_dz - T[i + 1, j] + T[i - 1, j] - (
(dr ** 2) / r[i]) * dT_dr + (
((-delta_H) * rho_catalyst * reaction_i * dr ** 2) / lambda_er))
# Boundary conditions
C_new[0, :] = C_new[1, :] # r=0 boundary
C_new[-1, :] = C_new[-2, :] # r=r_end boundary
T_new[0, :] = T_new[1, :] # r=0 boundary
T_new[-1, :] = Ts # r=r_end boundary
C_new[:, 0] = C0 # z=0 boundary
T_new[:, 0] = T0 # z=0 boundary
C_new[:, -1] = C_new[:, -2] # z=z_end boundary
T_new[:, -1] = T_new[:, -2] # z=z_end boundary
# convergence?
not_converged = np.max(np.abs(C_new - C)) > tolerance or np.max(np.abs(T_new - T)) > tolerance
# Calculate residuals
residual_C = np.max(np.abs(C_new - C))
residual_T = np.max(np.abs(T_new - T))
residuals_C.append(residual_C)
residuals_T.append(residual_T)
print(f"Residual for C in iteration {iteration}: {residual_C}")
print(f"Residual for T in iteration {iteration}: {residual_T}")
C = C_new
T = T_new
print(f"Iteration: {iteration + 1}")
print("Temperature:")
print(T)
print("Concentration:")
print(C)
print("------------------------")
iteration += 1
# create 2D grid for plot
R, Z = np.meshgrid(r, z, indexing='ij') # 修改为 'ij' indexing 以确保 R, Z 的形状与 T, C 匹配
plt.figure(figsize=(8, 6))
plt.subplot(121)
plt.contourf(Z, R, C, cmap='viridis')
plt.colorbar()
plt.title("Concentration Distribution")
plt.xlabel("z")
plt.ylabel("r")
plt.subplot(122)
plt.contourf(Z, R, T, cmap='viridis')
plt.colorbar()
plt.title("Temperature Distribution")
plt.xlabel("z")
plt.ylabel("r")
plt.tight_layout()
plt.show()
plt.figure(figsize=(10, 5))
plt.subplot(121)
plt.plot(residuals_C)
plt.yscale('log') # Use logarithmic scale for better visibility
plt.title("Residuals of C per Iteration")
plt.xlabel("Iteration")
plt.ylabel("Residual of C")
plt.subplot(122)
plt.plot(residuals_T)
plt.yscale('log') # Use logarithmic scale for better visibility
plt.title("Residuals of T per Iteration")
plt.xlabel("Iteration")
plt.ylabel("Residual of T")
plt.tight_layout()
