1D heat transfer

Python Help
1

Hi guys, I am in this forum and basically new at numerical modelling.
I am trying to solve a transient 1D heat transfer equation following a youtube tutorial and adapting it to my data. Even if there seem not to be any error, the code does not plot anything.
Can you help me to find the problem? Thank you a lot :slight_smile:

Here is my code:

import numpy as np
import matplotlib.pyplot as plt

L = 10
n = 10

T0= 0
T_rock = 600
T_magma = 1300

alpha = 10**(-6)

dx = L/n

t_final = 60
dt = 0.1

x = np.linspace(dx/2, L-dx/2, n)
t = np.arange(0, t_final, dt)

dTdt = np.empty(n)

T = np.ones(n)*T0

for j in range(1, len(t)):
for i in range(1, n-1):
dTdt[i] = alpha*(-(T[i]-T[i-1])/dx2 + (T[i+1]-T[i])/dx2)
dTdt[0] = alpha*(-(T0-T_magma)/dx2 + (T[1]-T0)/dx2)
dTdt[n-1] = alpha*(-(T[n-1]-T[n-2])/dx2 + (T_rock-T[n-1])/dx2)

T = T + dTdt*dt
plt.figure(1)
plt.plot(x,T)
plt.axis([0, L , 500, 1500])
plt.xlabel('Distance [m]')
plt.ylabel('Temperature [°C]')
plt.pause(0.05)

plt.show()

3

For a start I think you’ve missed declaring dx2 somewhere in your code. Check out my comment below:

import numpy as np
import matplotlib.pyplot as plt

L = 10
n = 10

T0= 0
T_rock = 600
T_magma = 1300

alpha = 10**(-6)

dx = L/n

t_final = 60
dt = 0.1

x = np.linspace(dx/2, L-dx/2, n)
t = np.arange(0, t_final, dt)

dTdt = np.empty(n)

T = np.ones(n)*T0

for j in range(1, len(t)):
    for i in range(1, n-1):
        dTdt[i] = alpha*(-(T[i]-T[i-1])/dx2 + (T[i+1]-T[i])/dx2)  # <---- where did dx2 come from?
        dTdt[0] = alpha*(-(T0-T_magma)/dx2 + (T[1]-T0)/dx2)
        dTdt[n-1] = alpha*(-(T[n-1]-T[n-2])/dx2 + (T_rock-T[n-1])/dx2)
        T = T + dTdt*dt

plt.figure(1)
plt.plot(x,T)
plt.axis([0, L , 500, 1500])
plt.xlabel('Distance [m]')
plt.ylabel('Temperature [°C]')
plt.pause(0.05)
plt.show()

Please post your traceback in tripple back tics for any subsequent errors.
They’re really helpful for us trying to help you :slight_smile:

2 Likes
4

I think in the original code it was dx**2. Markdown really messes the code up. I think that the forum should have some autodetection of (Python) code not between triple backticks. Almost all the new users here do not understand the problem with formatting.

5

I see it now. The ** followed by ** turned the mid section into bold instead.

6

There’s nothing wrong with the code. I’m running python3.9 on windows as this:

import numpy as np
import matplotlib.pyplot as plt

L = 10
n = 10

T0= 0
T_rock = 600
T_magma = 1300

alpha = 10**(-6)

dx = L/n

t_final = 60
dt = 0.1

x = np.linspace(dx/2, L-dx/2, n)
t = np.arange(0, t_final, dt)

dTdt = np.empty(n)

T = np.ones(n)*T0

for j in range(1, len(t)):
    for i in range(1, n-1):
        dTdt[i] = alpha*(-(T[i]-T[i-1])/dx**2 + (T[i+1]-T[i])/dx**2)
        dTdt[0] = alpha*(-(T0-T_magma)/dx**2 + (T[1]-T0)/dx**2)
        dTdt[n-1] = alpha*(-(T[n-1]-T[n-2])/dx**2 + (T_rock-T[n-1])/dx**2)
        T = T + dTdt*dt

the issue is in the plotting:

plt.figure()
plt.xlabel('Distance [m]')
plt.ylabel('Temperature [°C]')
# plt.axis([0, L , 500, 1500])  # <--- disabled.
plt.plot(x,T)
# plt.pause(0.05)  # <--- disabled.
plt.show()

where I get this after disabling two rows.

output

1 Like
7

Yes, that a was a naif error. I have found out where the problem was and, beyond the correction of dx**2, also parameter setting was lacking of physical meanng. Thank you very much for your support and suggestions! :slight_smile:

8

You’re welcome :wink:

1 Like