有限差分解2D PDE的矩阵初始化
2D稳态传热传质在圆柱坐标系中的PDE求解 用的是中心差分法
边界条件如下图所示
我应该怎样初始化我的网格？如果初始化是0的话 收敛到e-7后 边缘会是0 如下图所示
这里面有两个PDE 能量平衡（计算温度）和质量平衡（计算浓度）反应速率是温度和浓度的函数
我的代码里 我先用初始反应速率更新了每个网格的温度和浓度 然后更新了每个网格的反应速率再去得到新的温度和浓度
这样一直迭代温度和浓度
我这样的方法是对的吗？我矩阵初始化应该是都是0 还是都是初值呢？比如温度是570K 浓度是3622 mol/m3
以下是代码
import numpy as np
import matplotlib.ticker as ticker
import matplotlib.pyplot as plt
# Wall
# -----------
# | A D G |
# Inlet | B E H | Outlet
# | C F I |
# ------------
# Symmetric line
# Set total range
r_range = 0.01 # m
z_range = 0.25 # m
# Set grid size
grid_size = 0.0003 # 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.zeros((num_points_r, num_points_z))
# C = np.ones((num_points_r, num_points_z)) * 3622
T = np.zeros((num_points_r, num_points_z))
# T = np.ones((num_points_r, num_points_z)) * 570
Q_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.0000253856 # m2/s
rho_catalyst = 5000 # 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
Ts = 570 # wall temperature
A = u_s * rho_fluid * cp_fluid
# Set reaction parameter
k = 0.000005 # m3/kg s
E = 157000 # J/mol
R_gas = 8.314 # J/mol K
# Set initial parameter
C0 = 3622 # mol/m3 rho/MW 704 kg/m3 / 0.19435 kg/mol
T0 = 570 # K
# # 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[:, :] = k * np.exp(-E / (R_gas * T0)) * C0 # initial reaction rate
print("Initial reaction rate [1][0] ", Q_reaction[1, 1])
# iterative to convergence
tolerance = 1e-6
not_converged = True
iteration = 0
max_iterations = 10000
not_converged_2 = True
iteration_2 = 0
max_iterations_2 = 100
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)
C_new = np.copy(C)
Q_reaction_new = np.copy(Q_reaction)
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 / 4) * (T[i - 1, j] + T[i + 1, j] + T[i, j - 1] + T[i, j + 1]) + (dr / (8 * r[i])) * \
(T[i + 1, j] - T[i - 1, j]) - ((A * dr) / (8 * lambda_er)) * (T[i, j + 1] - T[i, j - 1]) + \
(-delta_H * 6 * Q_reaction_new[i, j] * rho_catalyst * dr ** 2) / (4 * lambda_er)
C_new[i, j] = (1 / 4) * (C[i - 1, j] + C[i + 1, j] + C[i, j + 1] + C[i, j - 1]) + (dr / (8 * r[i])) * (
C[i + 1, j] - C[i - 1, j]) - \
((u_s * dr) / (8 * Der)) * (C[i, j + 1] - C[i, j - 1]) - (
rho_catalyst * Q_reaction_new[i, j]) * \
dr ** 2 / (Der * 4)
not_converged = np.max(np.abs(T_new - T)) > tolerance or np.max(np.abs(C_new - C)) > tolerance
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)
C[0, :] = C[1, :] # Symmetric axis, lower boundary
C[-1, :] = C[-2, :] # Upper boundary
C[:, -1] = C[:, -2] # Outlet
T[0, :] = T[1, :] # Symmetric axis, lower boundary
T[:, -1] = T[:, -2] # Outlet
T = T_new
C = C_new
print(f"Residual for T in iteration {iteration}: {residual_T}")
print(f"Residual for C in iteration {iteration}: {residual_C}")
print(f"Iteration: {iteration + 1}")
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
Q_reaction_new[i, j] = k * np.exp(-E / (R_gas * T_new[i, j])) * C[i, j]
iteration += 1
fig, ax = plt.subplots(2, 1, figsize=(10, 10)) # 2 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: '{:.7f}'.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: '{:.7f}'.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()
plt.figure(figsize=(10, 6))
# Plot residuals_T
plt.plot(residuals_T, label='Residuals T')
# Plot residuals_C
plt.plot(residuals_C, label='Residuals C')
plt.title('Residuals')
plt.xlabel('Iteration')
plt.ylabel('Residual')
plt.legend()
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.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.ylabel("Concentration")
plt.legend()
plt.tight_layout()
plt.show()
- 以下回答由chatgpt基于相关博客总结生成:
问题:如何正确初始化用于解2D圆柱坐标系中稳态传热传质PDE的差分矩阵？
首先，让我们来了解一下问题的背景。您正在使用中心差分法来解决2D稳态传热传质问题，这是一个涉及温度和浓度的PDE（偏微分方程）。在求解过程中，您需要考虑能量平衡和质量平衡，并且反应速率是温度和浓度的函数。
在初始化差分矩阵时，您应该使用初始条件对矩阵进行初始化，而不是全部为0。将初始条件设置为0会导致在收敛到e-7后，边缘值变为0，这不符合实际情况。
更具体地说，您可以根据具体的问题和实际情况选择初始条件的值。例如，您可以使用温度为570K，浓度为3622mol/m3的值作为初始条件。这样可以更好地近似实际问题，并更快地收敛到正确的解决方案。
以下是一个示例代码，演示了如何初始化差分矩阵并使用初始条件：
def initialize_matrix(): rows = 10 # 矩阵的行数 cols = 10 # 矩阵的列数 # 创建一个空矩阵 matrix = [[0.0 for _ in range(cols)] for _ in range(rows)] # 根据初始条件更新矩阵的值 for i in range(rows): for j in range(cols): matrix[i][j] = initial_condition(i, j) # 使用初始条件函数获得初始值 return matrix def initial_condition(i, j): # 在这里根据具体问题的要求返回初始条件的值 # 例如，可以返回温度和浓度的初始值 temperature = 570 # 温度的初始值：570K concentration = 3622 # 浓度的初始值：3622mol/m^3 # 根据网格点的位置返回初始条件的值 # 例如，可以根据网格点的位置调整初始条件的值 return temperature + i + j, concentration - i - j # 调用初始化矩阵的函数 matrix = initialize_matrix() # 打印矩阵的值 for row in matrix: print(row)在上述示例代码中，
initialize_matrix函数创建了一个指定行数和列数的空矩阵，并使用initial_condition函数根据初始条件更新矩阵的值。initial_condition函数根据网格点的位置返回初始条件的值，您可以根据实际问题进行调整。请注意，以上示例代码只是一个演示，您需要根据具体问题进行修改和适应。
希望这个解决方案能够帮助您正确初始化用于解决2D圆柱坐标系中稳态传热传质PDE的差分矩阵！如果还有其他问题，请随时提问。
如果你已经解决了该问题, 非常希望你能够分享一下解决方案, 写成博客, 将相关链接放在评论区, 以帮助更多的人 ^-^