5
$\begingroup$

I would like to determine a temporal derivative and a Laplacian by the finite differences method to solve the heat equation in a 1-dimensional case. My aim is to get the sources term that is why I need to get all left terms of the equation.

When I plot the sources term as a function of time, the result is surprising... very far from expected. Is this a problem with the np.reshape? Thanks for help ;)

#profil1D_theta.shape[0] : time
#profil1D_theta.shape[1] : space

#profil1D_theta_reshape.shape                                                                                                                                                           
#Out[986]: (328, 411)

T_tempo = 1./115. #time period

#Time derivative

dt_theta_list = []

for i in range(profil1D_theta.shape[0]-1):
    for j in range(profil1D_theta.shape[1]-1):
        dt_theta = (profil1D_theta[i+1,j] - profil1D_theta[i,j])/T_tempo
        dt_theta_list.append(dt_theta)

dt_theta = np.asarray(dt_theta_list)
dt_theta_reshape = dt_theta.reshape(328,411)


#Laplacian

T_spatial = 0.000021661 #spatial period

dx2test_theta_list =[]
for i in range(0,profil1D_theta.shape[0]-1):
    for j in range(0,profil1D_theta.shape[1]-1):
        dxtest_theta = (profil1D_theta[i,j+1] - 2*profil1D_theta[i,j] + profil1D_theta[i,j-1])/T_spatial**2
        dx2test_list.append(dxtest_theta)

dx2test_theta = np.asarray(dx2test_theta_list)

dx2_theta_reshape = dx2test_theta.reshape(328,411)



rho = 7850. 
C = 500. 
k = 42. 
D = k/(rho*C) 
e = 3.*0.001 
b = 6.*0.001 
Boltz = 5.679373*10**(-8) 
Emissivite_echant = 0.97 #
T0 = 273.15 + 26. 
h = 5.

tau_th = rho*C*e/(2*h+8*Boltz*Emissivite_echant*T0**3)


#Sources

profil1D_theta_reshape = profil1D_theta[0:328,0:411]

w1_chaleur_list=[]

for i in range(0,profil1D_theta_reshape.shape[1]):
    for j in range(0,profil1D_theta_reshape.shape[0]):
        w_chaleur1D = rho*C*(dt_theta_reshape[j,i]+1./(tau_th)*profil1D_theta_reshape[j,i] - D*dx_theta_reshape[j,i])
        w1_chaleur_list.append(w_chaleur1D)

d,e = dt_theta_reshape.shape
w1_chaleur = np.asarray(w1_chaleur_list)
w1_chaleur_reshape = np.reshape(w1_chaleur,(e,d))

clim_min = 0
clim_max = 1*10**6

plt.imshow(w1_chaleur_reshape);plt.clim([clim_min,clim_max]);plt.colorbar();
plt.ylabel(r'Axe longitudinal')
plt.xlabel(r'Temps (images)')
plt.show()
Improve this question
$\endgroup$
5
  • 13
    $\begingroup$ You're not likely to get help with such a vague question and only hints of what's wrong. Please add more detail about your results (such as plots and what is "expected") and why you suspect the problem might be with np.reshape. As the question is written, you're basically asking us to just read and debug your code. $\endgroup$
    – Doug Lipinski
    Jul 17, 2014 at 23:34
  • $\begingroup$ You derivatives are edge-centered. This means that the size of the result array is going to be smaller than the input array. Also, you Laplacian code will likely do a bunch of things wrong given the lack of bounds checking. See numpy.diff() for a smarter way to approach this all. $\endgroup$
    – meawoppl
    Jul 20, 2014 at 23:58
  • $\begingroup$ I can add a bit of example code of how to do this if you are interested. $\endgroup$
    – meawoppl
    Jul 21, 2014 at 0:02
  • $\begingroup$ I'm coding my own version of Heat Equation Solver, in 2D. You might want to take a look at how I solved the Laplacian problem for the border of the grid. https://github.com/matpompili/HeatFlow $\endgroup$
    – Matteo Pompili
    Aug 15, 2014 at 22:04
  • $\begingroup$ You should check out my answer to this question about the same problem solving the heat equation. It's written for matlab but it should be easy to port to python. $\endgroup$
    – nluigi
    Mar 8, 2016 at 12:48

0

Reset to default

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Browse other questions tagged or ask your own question.