Reputation: 35
So I am trying to solve the heat equation using python.
The principle is that while the outdoor temperature varies through the day as so
T(0, t) = T0 + Tm cos(ω t) (boundary limit), well the temperature in the soil (4cm deep) also varies as a sinusoidal function, but the amplitude will be smaller (since the soil acts as insulation).
My code should plot the time on x and the temperature in the soil at 4cm deep on y.
My issue is that I do not have a sinusoidal function returned by Python
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 + 24 * 60 * 60 # End time in seconds (next day 6:50 am)
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))
# Set the initial conditions
temperature[:, 0] = 29.364 # Initial temperature at all depths (units: Celsius)
# Set the boundary condition at x = 0 (surface)
To = 29.364 # Initial temperature at 4 cm depth (units: Celsius)
Tm = 2.0 # Amplitude of external temperature variation (units: Celsius)
w = 2 * np.pi / (24 * 60 * 60) # Angular frequency of external temperature variation (units: rad/s)
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 the boundary condition at x = 0 (surface)
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
# Plotting
time = np.linspace(time_start, time_end, num_t_steps)
temperature_variation = temperature[int(depth / delta_x), :]
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')
plt.grid(True)
plt.show()
Upvotes: 1
Views: 83