提问人：Yovan Bach 提问时间：10/27/2023 最后编辑：sbottingotaYovan Bach 更新时间：10/29/2023 访问量：86
Python 中的一维热感应方程
1D Heat Induction Equation in Python
问：
我是 Python 的新手，我正在尝试求解 Python 中头部感应方程的解析解。我写了一个代码，但它没有按照我的意愿工作，我无法修复它。
在这里，您可以找到给定的解析解：
这是我的代码：
import numpy as np
import matplotlib.pyplot as plt
print("1D heat equation solver for a wall (Analytical Solution)")
wall_thickness = 31 # Thickness of the wall in centimeters
L = wall_thickness # Length of the wall
max_time = 1800 # Total simulation time in seconds
delta_x = 0.05
delta_t = 0.01
D = 93 #Diffusivity coeficient
# Define the spatial grid
x = np.arange(0, L + delta_x, delta_x)
# Initialize solution: the grid of u(x, t)
u = np.empty((max_time, len(x)))
# Set the boundary conditions
u_sur = 149 # Surrounding temperature
u_in = 38.0 # Initial temperature inside the wall
# Set the initial condition
u[0] = u_in
Tbase = u_sur + 2 * (u_in - u_sur)
def analytical_solution(x, t):
T = np.zeros_like(x)
for m in range(1, 100): # Sum from m=1 to a large number (e.g., 100)
T += np.exp(-D*((m * np.pi) / L)**2*t) * ((1 - (-1) ** m) / (m * np.pi)) * np.sin((m * np.pi * x) / L)
return Tbase + T
# Calculate the analytical solution for each time step
for k in range(1, max_time):
t = k * delta_t
u[k] = analytical_solution(x, t)
#print(f"Valeurs de la sol analytique {u[k]}.\n")
# Plot the temperature profiles at specific time points
plt.figure(figsize=(10, 6))
time_points = [0, 360, 720, 1080, 1440, 1800]
for t in time_points:
k = int(t / delta_t)
if k < max_time:
plt.plot(x, u[k], label=f"t = {t} seconds")
plt.title("Temperature Profiles (Analytical Solution) at Specific Time Points")
plt.xlabel("x (thickness)")
plt.ylabel("Temperature")
plt.ylim(0,150)
plt.xlim(0, L)
plt.legend()
plt.show()
print("Done!")
以下是我获得的结果：
正如你所看到的，结果是不连贯的，用红色绘制的曲线是我所期望的......
提前感谢您的帮助，
答：
0赞
lastchance
10/29/2023
#1
它是热传导，而不是热感应。
你有几个问题。首先，您似乎混淆了步骤数 （k） 和时间 （t）。
其次，正如 @slothrop 在评论中已经指出的那样，将总和与初始边界条件和边边界条件相结合时存在错误。
第三，看表达式 （1 - （-1） ** m）。如果 m 为奇数，则为 2，如果 m 为偶数，则为 0 - 因此，您可以通过仅考虑奇数 m 并将此项设置为 2 来大幅简化总和。
第四，一个物理问题 - 在给定时间内扩散的距离是sqrt（扩散率*时间）的数量级，你的扩散率为93个单位，并且必须扩散大约一半的墙宽度。如果这些“时间”（实际索引）点是你所说的，我建议你使用较小的时间步长。
试试这个：
import numpy as np
import matplotlib.pyplot as plt
print("1D heat equation solver for a wall (Analytical Solution)")
wall_thickness = 31.0 # Thickness of the wall in centimeters
L = wall_thickness # Length of the wall
max_step = 1800 # Number of timesteps
delta_x = 0.05
delta_t = 0.001
D = 93 # Diffusivity coeficient
# Define the spatial grid
x = np.arange(0, L + delta_x, delta_x )
# Initialize solution: the grid of u(x, t)
u = np.empty( ( 1 + max_step, len( x ) ) )
# Set the boundary conditions
u_sur = 149 # Surrounding temperature
u_in = 38.0 # Initial temperature inside the wall
# Set the initial condition
u[0] = u_in
def analytical_solution(x, t):
T = np.zeros_like(x)
for m in range( 1, 100, 2 ): # Sum from m=1 to a large number (e.g., 100)
T += np.exp( - D * ( m * np.pi / L ) ** 2 * t ) * 2.0 / ( m * np.pi ) * np.sin( m * np.pi * x / L )
return u_sur + 2 * ( u_in - u_sur ) * T
# Calculate the analytical solution for each time step
for k in range( 1, 1 + max_step ):
t = k * delta_t
u[k] = analytical_solution( x, t )
#print(f"Valeurs de la sol analytique {u[k]}.\n")
# Plot the temperature profiles at specific time points
plt.figure(figsize=(10, 6))
time_points = [ 0, 360, 720, 1080, 1440, 1800 ]
for k in time_points:
if k <= max_step:
plt.plot(x, u[k], label=f"t = { k * delta_t } seconds")
plt.title("Temperature Profiles (Analytical Solution) at Specific Time Points")
plt.xlabel("x (thickness)")
plt.ylabel("Temperature")
plt.ylim(0,150)
plt.xlim(0, L)
plt.legend()
plt.show()
print("Done!")
输出：
评论
0赞
Yovan Bach
11/5/2023
非常感谢您的帮助和这些很好的解释！
下一个：追加以列出 JAX 的替代工作流
评论
Σ ...2(T_in - T_sur)T_sur + 2(T_in - T_sur) + Σ...