急 Python二维稳态传热传质 偏微分方程
求指导 可以加钱！
固定床催化反应，二维稳态传热传质问题，有源项，解线性偏微分方程
请问我这个离散做得对吗
传热方程 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()
采用chatgpt：
您的代码看起来基本上是正确的。我注意到在计算反应速率时，您使用了温度T[i, j]，而不是新的温度T_new[i, j]。为了正确计算反应速率，您需要使用更新后的温度值T_new。所以您需要将以下行：
reaction_i = k * np.exp(-E / (R_gas * T[i, j])) # m3/kg s
改为：
reaction_i = k * np.exp(-E / (R_gas * T_new[i, j])) # m3/kg s
另外，我注意到在计算C_new和T_new时，您的代码中的一些项位置不正确。在更新C_new时，您需要将rho_catalyst * reaction_i乘以r[i]，而不是除以r[i]。此外，您在更新T_new时的括号位置也有些问题。请参考以下更正的代码段：
# 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 * r[i] * dr ** 2) / lambda_er)
这些更正应该能帮助您解决收敛问题。请注意，如果问题仍然存在，您可能需要调整迭代参数（如tolerance和max_iterations）或初始条件（如C0和T0）以更好地满足收敛要求。
可以帮我修改代码吗 或者用别的方法来实现
你没问ChatGPT吗
运行了你的代码了，只有一些警告
reaction_i = k * np.exp(-E / (R_gas * T[i, j])) # m3/kg s
改为：
reaction_i = k * np.exp(-E / (R_gas * T_new[i, j])) # m3/kg s
试一下
Python二维稳态传热传质 偏微分方程
非常详细，可以参考下
可能是的离散化步长的选择问题，如果离散化步长过大，会导致数值误差较大，影响数值解的准确性。如果离散化步长过小，会导致计算量过大，计算时间过长。
或者选择了不合适的边界条件。中心差分法是一种数值方法，但是它的精度可能不够高，可能会导致数值解不准确。可以尝试其他数值方法，例如有限元法或有限差分法。
检查一下初始条件
增加迭代次数试试
在固定床催化反应中，对二维稳态传热传质问题进行离散化时，遇到收敛困难，常见问题可能源于网格分辨率不够细致、选择不合适的迭代算法和处理边界条件错误，优化建议包括增加网格密度、尝试不同的迭代算法以及检查边界条件的正确性
引用chatgpt回答: 在进行二维稳态传热传质问题的解线性偏微分方程时遇到了问题。注意到您将中心差分方法应用于径向和轴向的偏微分方程，并且进行了迭代以达到收敛。以下是对您的代码的一些观察和建议：
反应速率计算问题：
- 在每次迭代中，您计算了反应速率reaction_i，但似乎没有考虑到温度变化（T）在反应速率中。您只使用了固定的温度（Ts）来计算reaction_i。您需要将当前位置（i，j）处的温度（T）用于每次迭代中的反应速率计算。
界面条件问题：
- 在更新C_new和T_new之前，您在边界条件上设置了C_new和T_new的值。然而，在边界条件的处理方面可能出现了一些问题。在r = 0边界处，应该选择适当的边界条件，而不是简单地将C_new和T_new的第一行或最后一行设置为相邻位置的值。
循环范围问题：
- 在更新C_new和T_new的循环中，您的r和z的循环从1开始，而不是从0开始。这可能会导致计算中缺少一些数据点。确保循环从0开始，从而包括所有的数据点。
边界条件处理和收敛条件：
- 您需要仔细检查边界条件的设置是否正确，并确保它们与问题的物理特性相匹配。另外，请检查收敛条件的设置是否合适，以便适应问题的求解。
根据上述观察和建议，您可以尝试对代码进行修改并重新运行，以看是否能够解决问题。同时，建议您仔细检查方程、边界条件和初始条件以确保符合问题的要求。如果问题仍然存在，请提供更详细的错误描述，以便我能够为您提供更具体的帮助。
我看到您在解决固定床催化反应的二维稳态传热传质问题时使用了中心差分法来离散化偏微分方程。从代码上看，您的离散化和迭代过程似乎是正确的，但可能出现了一些问题导致无法收敛。
以下是一些可能导致问题的原因：
迭代次数不够：您设置了最大迭代次数为1000次，但如果残差未能在这个迭代次数内收敛，可能需要增加迭代次数。
初始条件选择：初始条件对于迭代结果的收敛至关重要。请确保初始条件满足问题的物理约束，并尽量接近最终稳态的解。
离散化步长选择：检查离散化步长的选择是否合适。如果步长选择不当，可能会导致数值不稳定或者无法收敛。可以尝试减小步长并重新运行代码。
边界条件处理：确保边界条件的处理正确无误。特别是在计算梯度时，边界处的处理非常重要。
反应速率项：检查反应速率项的计算是否正确，包括温度项、气体常数等参数是否按照正确的单位和公式计算。
你的离散化过程看起来基本正确，但我注意到一些潜在的问题，可能导致你的收敛性问题：
- 传热方程的边界条件：在更新温度的循环中，你设置了 T_new[-1, :] = Ts，即 r=r_end 处的温度为固定的 Ts。然而，在边界条件中，你应该考虑热传导过程中的对流热传递。因此，你可能需要根据对流换热系数 Ut 和壁面温度 Ts 来更新边界条件。
- 更新温度的计算式：在更新温度的计算式中，有一项 (1/2) 的系数，并且有一个额外的加法项 ((-delta_H) * rho_catalyst * reaction_i * dr ** 2) / lambda_er。请确保这些系数和项的来源和推导正确，并与你的物理模型相匹配。
此外，你可以尝试以下方法来改进收敛性：
- 检查初始条件：确保初始的浓度和温度场满足问题的物理要求，并且它们在适当的范围内。
- 调整迭代次数和收敛准则：增加最大迭代次数，并适当调整收敛准则，以允许更多的迭代过程来达到收敛。
- 调整离散化步长：尝试减小 dr 和 dz 的值，以增加空间网格的分辨率，并检查结果是否收敛。
- 调整参数和常数：对于反应速率常数 k、活化能 E 和其他物理参数，确保它们的值正确并与问题的实际情况相匹配。
最后，确保对你的模型进行仔细验证，并与问题的实际情况进行比较和分析。如果问题仍然存在，可能需要进一步调试和排除故障，例如使用调试输出或绘制其他关键变量的变化情况，以便更好地理解问题所在。