# -*- coding: utf-8 -*-
"""
Created on Mon Nov 20 20:37:28 2023
@author: THOMAS-SCBoduch
"""
# plot a nice polar function
# Try with animation
import math
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.animation import FuncAnimation
num_r=120
num_theta=25
max_iter=2000
radius=1
delta_r = radius / num_r
delta_theta = 2*np.pi / num_theta
alpha=1.1
delta_t = (delta_r ** 2)/(4 * alpha)
gamma = (alpha * delta_t) / (delta_r ** 2)
##########################Def of the grid#############################
r=np.linspace(0, radius, num_r)
theta=np.linspace(0,2*np.pi, num_theta)
R,Theta = np.meshgrid(r,theta)
#T=np.ones((num_theta,num_r))
T=np.ones((num_theta,num_r))
######################################################################
####Def of the plot
fig, ax1 = plt.subplots(subplot_kw=dict(projection='polar'))
#plt.rcParams['pcolor.shading'] ='nearest'
quad = ax1.pcolormesh(Theta,R,T,shading='gouraud')
######################################################################
#other setups,
#ax1.plot(Theta, R , T, color='r', ls='none', marker='.')
#def des fonctions initiales
#conditions aux limites
T[:,-1]=20.0
T[:,0]=25.0
T[:,1]=23.0
#T[0,0]=25.0
#T[0,1]=23.0
#T[-1,0]=25.0
#T[-1,1]=23.0
#Conditions initiales
def animate(frame):
global T, quad
T_temp = T.copy() # Initialisez T_temp en dehors de la boucle
for i in range(1, num_r - 1):
r_i = i * delta_r
T_temp[:, i] = T[:, i] + delta_t * ((alpha / delta_r**2) * (T[:, i + 1] - 2 * T[:, i] + T[:, i - 1]) + (alpha / r_i / delta_r) * (T[:, i + 1] - T[:, i]))
#print(T_temp[:,1])
# Copiez les conditions initiales de la frame précédente
T_temp[:, 0] = T[:, 0]
T_temp[:, 1] = T[:, 1]
T_temp[:, -1] = T[:, -1]
#T_temp[0, 0] = T[0,0]
#T_temp[0, 1] = T[0,1]
#T_temp[-1, 0] = T[-1,0]
#T_temp[-1, 1] = T[-1,1]
T = T_temp.copy()
#quad.set_array(T.flatten())
# Mise à jour du plot
ax1.clear() # Effacer l'axe pour éviter la superposition des images
quad = ax1.pcolormesh(Theta, R, T, shading='gouraud', cmap=plt.cm.jet, vmin=-10 , vmax=50)
ax1.set_title('{:.1f} ms'.format(frame*delta_t*1000))
return quad,
""" Laplacien steady state R
#champ scalaire de temperature
T=np.ones((num_theta,num_r))
#conditions aux limites
T[:,-1]=20.0
#Conditions initiales
T[:,0]=25.0
T[:,1]=23.0
for i in range(1,num_r-1):
r_i = i * delta_r
T[:,i+1]= (T[:,i] * (delta_r + 2*r_i) - r_i * T[:,i-1] ) / ( r_i + delta_r )
print(T[0,i+1])
"""
""" Laplacien steady state R and Theta
#champ scalaire de temperature
T=np.ones((num_theta,num_r))
#conditions aux limites
T[:,-1]=20.0
#Conditions initiales
T[:,0]=25.0
T[:,1]=23.0
for i in range(1,num_r-1):
r_i = i * delta_r
for j in range(1,num_theta-1):
T[:,i+1]= (T[:,i] * (delta_r + 2*r_i) - r_i * T[:,i-1] ) / ( r_i + delta_r )
T[j+1,:]= 2*T[j,:] - T[j-1,:]
#"""
anim=FuncAnimation(fig,animate,interval=50, frames=max_iter,repeat=False)
anim.save("heat_equation_solution.gif")
plt.colorbar(quad)
plt.show()
Initial condition are as follow :
The center of the cylinder is at 25°C (this is where the gas filling is going on). The surface of the cylinder is at 20°C (ambiant temperature)