Python - 通过有限差分计算时间导数和laplacien

计算科学 有限差分 Python 麻木的
2021-12-23 08:19:22

我想通过有限差分法确定时间导数和拉普拉斯算子,以解决一维情况下的热方程。我的目标是得到源项,这就是为什么我需要得到方程的所有左项。

当我将源项绘制为时间的函数时,结果令人惊讶……与预期相去甚远。这是 np.reshape 的问题吗?感谢帮助 ;)

#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()
0个回答
没有发现任何回复~