Elisa909
Elisa909

Reputation: 35

Heat equation plot issues

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

Answers (0)

Related Questions