圆柱坐标二维稳态控制偏微分方程离散化
我正在开发一个固定床反应器中的热量和质量传递模型，这个反应器中存在吸热反应。
如下所示的是一个二维图，从左侧有恒定温度和浓度的入口，顶部有恒定温度的墙壁，底部是一个中心线（对称线），温度和浓度没有变化，出口（右侧）也一样，这些都定义了初始条件和边界条件。
存在三个方程，即质量平衡，能量平衡和反应速率，两个控制方程是通过反应速率耦合在一起的。
我计划使用Python中的有限差分方法来解这些偏微分方程。首先，我使用中心差分近似法来离散化这些方程。
质量平衡：
能量平衡：
反应速率没有离散化，因为它使用的是当前的T（温度）和C（浓度）。
问题1：当我离散化能量平衡方程时，我需要将所有的T[i,j]移至左侧作为未知数吗？因为能量平衡包括了一个反应速率项，它是T和C的函数，反应项中的T就是T[i,j]。我现在的做法是，将反应项留在右侧，因为这个温度是节点中前一个的温度。这里我有些困惑。
问题2：我从使用初始条件计算T开始循环，循环将按列逐列开始。例如，首先将使用其周围四个节点来计算T[1, 1]，其中有两个节点是已知值，然后将计算T[2,1]，直到T[3,1]。底部和出口的边界条件是没有变化，所以我设定T[-1,:] = T[-2,:]作为底部的条件，和T[:,-1] = T[:,-2]作为出口的条件。然而，由于数组是以0值初始化的，这些边界条件会影响结果，当靠近边界时，值会减小至0。那么，我应该如何在Python中正确设置Neumann BC（边界条件）呢？
问题3：你能给一些关于如何解这三个方程的建议吗？目前，我有一个循环并按序解这三个方程。我用初始值得到一个新的温度节点，然后使用这个温度和初始浓度得到该节点的反应速率。新的反应速率用于得到更新的浓度。关于这个迭代，我感到有些困惑。
我已经在下面附上了我的代码，不过结果并不正确。谢谢你的建议和指导！
import numpy as np
import matplotlib.pyplot as plt
# Set total range
r_range = 10 # mm
z_range = 250
# Set grid size
grid_size = 1
# 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.zeros((num_points_r, num_points_z))
T = np.zeros((num_points_r, num_points_z))
# T = np.ones((num_points_r, num_points_z)) * 570
reaction = 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
Der = 0.001
rho_catalyst = 4000 # kg/m3
rho_fluid = 704 # kg/m3
cp_fluid = 3037 # J/ kg K
lambda_er = 39.8 # 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
A = u_s * rho_fluid * cp_fluid / 1000
# Set reaction parameter
k = 0.0108 # m3/kg s
# k = 0.00008 # m3/kg s
E = 117000 # J/mol
R_gas = 8.314 # J/mol K
H = 117000 # J/ mol H2
# Set initial parameter
C0 = 5134 # mol/m3 rho/MW kg/m3 / kg/mol
T0 = 570 # K
# # Boundary conditions
T[:, 0] = T0 # inlet temp, left 570K
T[:, -1] = 25 # outlet right
T[-1, :] = Ts # constant wall temp, up boundary
T[0, :] = 300 # lower boundary
C[:, 0] = C0 # inlet concentration
C[-1, :] = T[-2, :]
C[0, :] = C[1, :]
reaction[:, 0] = reaction[-1, :] = k * np.exp(-E / (R_gas * T0)) * C0 # initial reaction rate
print("Initial reaction rate ", reaction)
# iterative to convergence
tolerance = 1e-7
not_converged = True
iteration = 0
max_iterations = 1000
print("A:", A)
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 new temperature
T_new[i, j] = (1 / 2) * (T[i - 1, j] + T[i + 1, j]) + (dr / (4 * r[i])) * (T[i + 1, j] - T[i - 1, j]) - \
((A * dr) / (4 * lambda_er)) * (T[i, j + 1] - T[i, j - 1]) - \
(H * reaction[i, j] * rho_catalyst * dr ** 2) / (2 * lambda_er)
# T_new[i, j] = (1 / (8 * r[i])) * (dr * (T[i + 1, j] - T[i - 1, j]) + 2 * r[i] * (
# T[i - 1, j] + T[i + 1, j] + T[i, j - 1] + T[i, j + 1])) - (A / (lambda_er * 8 * r[i])) * (
# dr * (T[i + 1, j] - T[i - 1, j])) - (H * reaction * rho_catalyst) / (8 * A * r[i])
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 using average of old and new temperature
reaction[i, j] = k * np.exp(-E / (R_gas * (T_new[i, j - 1] + 1e-10))) * C[i, j - 1]
# Calculate new concentration
C_new[i, j] = (1 / 2) * (C[i - 1, j] + C[i + 1, j]) + (dr / (4 * r[i])) * (C[i + 1, j] - C[i - 1, j]) - \
((u_s * dr) / (4 * Der)) * (C[i, j + 1] - C[i, j - 1]) + (rho_catalyst * reaction[i, j]) / (
Der * 4 * r[i])
# Update temperature with reaction heat
# T_new[i, j] -= (H * reaction * rho_catalyst) / (4 * r[i] * lambda_er)
# Check for convergence
not_converged = np.max(np.abs(T_new - T)) > tolerance or np.max(np.abs(C_new - C)) > tolerance
# Calculate residuals
residual_T = np.max(np.abs(T_new - T))
# residual_C = np.max(np.abs(C_new - C))
residuals_T.append(residual_T)
# residuals_C.append(residual_C)
print(f"Residual for T in iteration {iteration}: {residual_T}")
# print(f"Residual for C in iteration {iteration}: {residual_C}")
# Update C and 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.contourf(Z, R, T, levels=100, cmap='viridis')
plt.colorbar()
plt.title("Temperature Distribution")
plt.xlabel("z")
plt.ylabel("r")
plt.figure(figsize=(8, 6))
plt.contourf(Z, R, C, levels=100, cmap='viridis')
plt.colorbar()
plt.title("Concentration Distribution")
plt.xlabel("z")
plt.ylabel("r")
plt.tight_layout()
plt.show()
plt.figure(figsize=(10, 5))
plt.plot(residuals_T)
plt.title("Residuals of T per Iteration")
plt.xlabel("Iteration")
plt.ylabel("Residual of T")
plt.tight_layout()
plt.show()
z_index = num_points_z // 2 # Index corresponding to z = z_range/2
plt.figure(figsize=(10, 5))
# Subplot for z=0
plt.subplot(1, 3, 1) # 1 row, 2 columns, first subplot
plt.plot(r, T[:, 0], label=f'z = {z_start}')
plt.title("Temperature Trend at z = 0")
plt.xlabel("r")
plt.ylabel("Temperature")
plt.legend()
# Subplot for z=z/2
plt.subplot(1, 3, 2) # 1 row, 2 columns, first subplot
plt.plot(r, T[:, z_index], label=f'z = {z_start}')
plt.title("Temperature Trend at z = z_range/2")
plt.xlabel("r")
plt.ylabel("Temperature")
plt.legend()
# Subplot for z=250
plt.subplot(1, 3, 3) # 1 row, 3 columns, second subplot
plt.plot(r, T[:, -1], label=f'z = {z_end}')
plt.title("Temperature Trend at z = 250")
plt.xlabel("r")
plt.ylabel("Temperature")
plt.legend()
plt.tight_layout()
plt.show()
- 这个问题的回答你可以参考下: https://ask.csdn.net/questions/254897
圆柱体稳态温度分布 | 分离变量法（二）| 偏微分方程（十四）
可以参考这个例子
https://blog.51cto.com/u_13341/6590004
具体报什么错误了么，把错误提示贴上来，根据报错信息才好找原因
可以将所有的T[i,j]移至左侧作为未知数，这样可以将方程转化为求解T[i,j]的形式。反应速率项可以是T[i,j]和C的函数，因此可以将其保留在右侧。确保你的方程在左侧和右侧的项之间平衡，并满足边界条件。
在Python中，可以使用索引来设置边界条件。
合理的设置边界条件