热方程绘图异常：Python求解土壤温度未生成正弦曲线

问题分析 #

你的代码存在两个核心问题导致结果不符合预期：

• 最深处边界条件未处理：代码仅更新了中间层（x从1到num_x_steps-2）的温度，最深处（x=4cm）的温度始终保持初始值，完全没有参与热传导过程，破坏了温度场的连续性。
• 初始瞬态干扰：初始时刻所有深度温度统一设为29.364，但表面边界条件在t=0时的温度是To + Tm*cos(w*time_start)，两者大概率不一致，系统需要时间过渡到稳态正弦响应，仅模拟24小时不足以让瞬态完全消失。

修复后的代码 #

import numpy as np
import matplotlib.pyplot as plt

# Constants
D = 4e-7  # Thermal diffusivity of the soil (units: cm^2/s)
depth = 4.0  # Depth in cm
time_start = 6 * 60 * 60 + 50 * 60  # Start time in seconds (6:50 am)
time_end = time_start + 72 * 60 * 60  # 延长模拟时间至72小时，确保系统达到稳态
delta_t = 600  # Time step size in seconds (10 minutes)

# Discretization
delta_x = 0.1  # Spatial step size in cm

# Calculate the number of spatial and temporal steps
num_x_steps = int(depth / delta_x) + 1
num_t_steps = int((time_end - time_start) / delta_t) + 1

# Initialize temperature array
temperature = np.zeros((num_x_steps, num_t_steps))

# 用热方程稳态解析解初始化，消除初始瞬态干扰
To = 29.364  # Mean temperature
Tm = 2.0  # Amplitude of external temperature variation (units: Celsius)
w = 2 * np.pi / (24 * 60 * 60)  # Angular frequency
sqrt_term = np.sqrt(w / (2 * D))
x_vals = np.linspace(0, depth, num_x_steps)
initial_t = time_start
temperature[:, 0] = To + Tm * np.exp(-x_vals * sqrt_term) * np.cos(w * initial_t - x_vals * sqrt_term)

# Surface boundary condition
boundary_condition = To + Tm * np.cos(w * np.linspace(time_start, time_end, num_t_steps))

# Solve the heat equation
for t in range(1, num_t_steps):
    temperature[0, t] = boundary_condition[t]  # Update surface temperature
    # Update middle layers
    for x in range(1, num_x_steps - 1):
        temperature[x, t] = temperature[x, t - 1] + D * (
            (temperature[x + 1, t - 1] - 2 * temperature[x, t - 1] + temperature[x - 1, t - 1]) / (delta_x ** 2)
        ) * delta_t
    # 处理最深处绝热边界（∂T/∂x = 0）
    temperature[-1, t] = temperature[-2, t]

# 仅绘制最后24小时的稳态数据
plot_start_idx = num_t_steps - int(24*60*60/delta_t)
time = np.linspace(time_start + (time_end - time_start)*(plot_start_idx/num_t_steps), time_end, num_t_steps - plot_start_idx)
temperature_variation = temperature[int(depth / delta_x), plot_start_idx:]

plt.plot(time, temperature_variation, linestyle='-', color='blue')
plt.xlabel('Time (s)')
plt.ylabel('Temperature (°C)')
plt.title('Temperature Variation in Soil at 4 cm Depth (Steady-State)')
plt.grid(True)
plt.show()

关键修复说明 #

• 补充最深处边界条件：设置土壤最深处为绝热边界（深层温度变化可忽略），通过temperature[-1, t] = temperature[-2, t]实现，保证热传导的连续性。
• 优化初始条件：使用热方程的解析稳态解初始化温度场，直接让系统接近稳态，大幅减少初始瞬态的影响。
• 延长模拟时间并筛选数据：模拟72小时后仅取最后24小时的数据绘制，确保结果是稳态的正弦响应——此时4cm深度的温度振幅会小于表面，且存在明显相位滞后，完全符合土壤隔热的物理规律。

内容的提问来源于stack exchange，提问作者Elisa909