Python 中的一维热感应方程

Question

我是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!")

以下是我得到的结果：

正如你所看到的结果是不连贯的，红色绘制的曲线正是我所期望的......

预先感谢您的帮助，

Answer 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!")

输出：