Il faudrait savoir dans quel cadre la question se pose, parce au vu de la géométrie l'écriture dun code à partir de zéro parait représenter une tache considérable.
A mon avis (mais c'est uniquement le mien...), cela ressort de l'utilisation d'un logiciel de type COMSOl ANSYS SimScale ...
On pourrait après comparer avec l'approximation Python.
J'ai essayé de faire quelque chose comme ça, rien ne se passe la complexité est désatreuse. Auriez-vous une solution ?
Code:image_disque = Image.open(r"C:\Users\Téo\Desktop\tipe disque.png") ar = array(image_disque) # Définition de la matrice: duree_simu = 1000 v = np.empty((duree_simu, np.shape(image_disque)[0], np.shape(image_disque)[1])) # conditions initiales: T_ini2 = 46.8 + 273.15 #K T_air2 = 25.3 + 273.15 #K T_metal = 35 + 273.15 #K #milieu du cercle: x_mi = np.shape(image_disque)[0]/2 y_mi = np.shape(image_disque)[1]/2 for i in range(0,np.shape(image_disque)[0]): for j in range(0,np.shape(image_disque)[1]): pixel = image_disque.getpixel((j, i)) if pixel == (0, 247, 64, 255): v[0,i,j] = T_air2 elif np.sqrt((i-x_mi)**2+(j-y_mi)**2)>=400 and np.sqrt((i-x_mi)**2+(j-y_mi)**2)<=448 and pixel != (0, 247, 64, 255) : v[0,i,j] = T_ini2 else: v[0,i,j] = T_metal def calculate(v): for k in range(0, duree_simu - 1): for r in range(1, round(np.shape(image_disque)[0]/2) ): for i in range(0,np.shape(image_disque)[0]-1): for j in range(0,np.shape(image_disque)[1]-1): if round(np.sqrt((i-x_mi)**2+(j-y_mi)**2)) == r : v[k + 1, i , j] = (dt*lambd / ( m*r*rho * C_th)) * (( (v[k, i,j ] - v[k, i,j]) / (2 * dr)) + r * m * (v[k, i,j] + v[k, i,j] - 2 * v[k, i,j]) / ( dr**2)) -(dt/(e * rho * C_th)) * (h * (v[k, i,j] - T_air) + epsilon * sigma * (v[k, i,j] ** 4 - T_air ** 4)) +v[k,i,j] return (v) v = calculate(v) plt.imshow(v[0], cmap='viridis', origin='lower', extent=(0, np.shape(image_disque)[1], 0, np.shape(image_disque)[0])) # Définissez les valeurs personnalisées pour l'axe des ordonnées et abscisses y_ticks = np.linspace(0, np.shape(image_disque)[0], 9) # Crée 9 valeurs de 0 à la hauteur de l'image x_ticks = np.linspace(0, np.shape(image_disque)[1], 9) # Crée 9 valeurs de 0 à la largeur de l'image plt.yticks(y_ticks, [-8, -6, -4, -2, 0, 2, 4, 6, 8]) # Personnalisez les étiquettes de l'axe des ordonnées plt.xticks(x_ticks, [-8, -6, -4, -2, 0, 2, 4, 6, '8 cm']) # Personnalisez les étiquettes de l'axe des abscisses plt.colorbar(label='Température (K)') plt.show()
Pardon, j'avais pas vu votre messsage. On ne pourrait donc pas adapter le code précédent ?
On pourrait l'adapter mais cela deviendrait une usine à gaz.
Pour ce qui est de votre code, récupérez la forme et initialiser à partir d'une photo est une bonne idée.
Mais il faut l'utiliser et pas seulement pour initialiser la température.
Il ne faut appliquer la loi de Fourier que là il y a de la matière.
Dans calculate, il ne devrait y avoir que trois boucles k i j : la boucle en r ne sert à rien (et le test correspondant) : vous n'avez que deux dimensions d'espace.
Par contre, il devrait y avoir un test pour savoir si on est bien sur le métal.
Vous pouvez continuer sur cette idée, cela reviendra à imposer comme CL T=Tair sur les faces perpendiculaires au dessin, très douteux mais simple.
Remarque : une équation avec une ligne qui déborde de l'écran rend les choses illisibles ; faites des calculs intermédiaires : flux convectif / radiatif les parmètres à regrouper et calculer une fois pour toutes...
Dit en clair et écrit en TEX (iliisible en Python)
On pose (et calcule une fois pour toutes) :
(on calcule en cartésiennes)
Dans la boucle :
Aux erreurs de calcul près ...
Je vous remercie pour votre réponse. J'ai modifié mon code, pour l'accélérer j'ai décide de diminuer la résolution de mon image. J'ai également fait quelque chose pour ne pas modifier la température des pixels appartenant au thermostat. Je crains cependant que cela impacte grandement la complexité du programme.
J'ai également un problème au niveau de la modélisation, à k = 5 les températures augmentent à 1200K, je ne comprends pas pourquoi. L'image utilisée si jamais vous souhaitiez tester de votre côté tipe disque.png
Figure 2023-10-26 135032.png
Dites moi ce que vous en pensez
Code:### Deuxième version # On compresse l'image pour optimiser le script image_disque_ini = Image.open(r"C:\Users\Téo\Desktop\tipe disque.png") dim_ini = image_disque_ini.size red = 5 # facteur de reduction image_disque =image_disque_ini.resize((int(dim_ini[0]/red), int(dim_ini[1]/red))) ar = array(image_disque) # Définition de la matrice: duree_simu = 10 v = np.empty((duree_simu, np.shape(image_disque)[0], np.shape(image_disque)[1])) # conditions initiales: T_ini2 = 46.8 + 273.15 #K T_air2 = 25.3 + 273.15 #K T_metal = 35 + 273.15 #K #milieu du cercle: x_mi = np.shape(image_disque)[0]/2 y_mi = np.shape(image_disque)[1]/2 vide = [] for i in range(0,np.shape(image_disque)[0]): for j in range(0,np.shape(image_disque)[1]): pixel = image_disque.getpixel((j, i)) if pixel == (0, 247, 64, 255): v[0,i,j] = T_air2 vide.append((i,j)) elif np.sqrt((i-x_mi)**2+(j-y_mi)**2)>=400/red and np.sqrt((i-x_mi)**2+(j-y_mi)**2)<=448/red and pixel != (0, 247, 64, 255) : v[0,i,j] = T_ini2 else: v[0,i,j] = T_air2 dx = 0.1 A = (dt*lambd / ( rho * C_th * m * dx**2)) B = (2*dt*h) / (e*rho*C_th) C = (2*dt*epsilon * sigma ) / (e*rho*C_th) def calculate(v): for k in range(0, duree_simu - 1): for i in range(0, np.shape(image_disque)[0]-1): for j in range(0, np.shape(image_disque)[1]-1): if (i,j) not in vide: v[k+1, i, j] = v[k, i, j] + A * (v[k, i+1, j] + v[k, i-1, j] + v[k, i, j+1] + v[k, i, j-1]) / 4 - (B * (v[k, i, j] - T_air2) + C * (v[k, i, j]**4 - T_air2**4)) return v v = calculate(v)
Pour l'affichage :
Code:plt.imshow(v[5], cmap='viridis', origin='lower', extent=(0, np.shape(image_disque)[1], 0, np.shape(image_disque)[0])) # Définissez les valeurs personnalisées pour l'axe des ordonnées et abscisses y_ticks = np.linspace(0, np.shape(image_disque)[0], 9) # Crée 9 valeurs de 0 à la hauteur de l'image x_ticks = np.linspace(0, np.shape(image_disque)[1], 9) # Crée 9 valeurs de 0 à la largeur de l'image plt.yticks(y_ticks, [-8, -6, -4, -2, 0, 2, 4, 6, 8]) # Personnalisez les étiquettes de l'axe des ordonnées plt.xticks(x_ticks, [-8, -6, -4, -2, 0, 2, 4, 6, '8 cm']) # Personnalisez les étiquettes de l'axe des abscisses plt.colorbar(label='Température (K)') plt.show()
Ok j'ai reperé le problème, j'ai introduit un facteur un m dans "A = (dt*lambd / ( rho * C_th * m * dx**2))". Après supression j'ai. Toutefois,ça me parrait étrange que tout mon disque soit à la même température en si peu de temps ( sans compter la bande chauffer par frottements), à noter que j'ai mis le disque (hors parties chauffée) à la température de l'air à t=0 .
Quelques remarques:
if pixel == (0, 247, 64, 255): dangereux à moins d'être sûr de votre fichier le 64 pourrait facilement être 63 ou 65
duree_simu = 10 est plutôt nbre_simu et duree_simu = nbre_simu*dt
dx=0.1 : vous ne pouvez pas imposé dx, c'est le pas de la grille dx=(largeur réelle (16 cm ?))/(largeur en pixel)
Pour alléger, j'ai calculé
nX= np.shape(image_disque)[0]
nY = np.shape(image_disque)[1]
puis remplacer dans la suite
for i in range(0, nX-1): OK pour nX-1 mais il faut faire la même chose au début donc remplacer 0 par 1 sinon v[k, i-1, j] va poser pb en i=0
Le temps caractéristique de convection est 300 s, de conduction 1300 s ; si cela va nettement plus vite, c'est qu'il y a un pb :
"imposer comme CL T=Tair sur les faces perpendiculaires au dessin, très douteux mais simple" est donc tellement douteux que cela ne marche pas.
Il faudrait donc appliquer des CL raisonnable sur ces faces mais cela deviendrait vraiment usine à gaz.
Très bien j'ai modifié, pour accélérer le code ma liste "vide" n'est plus une liste python mais une liste numpy, c'est nettement plus rapide.
Code:J'utilise également la fonction njit de numba, je ne sais pas trop ce que ça vaut mais ça à l'air OK. Cependant les températures s'emballent, 8000 K !!!. Je ne comprends pas trop pourquoi ### Deuxième version # On compresse l'image pour optimiser le script image_disque_ini = Image.open(r"C:\Users\Téo\Desktop\tipe disque.png") dim_ini = image_disque_ini.size red = 3 # facteur de reduction image_disque =image_disque_ini.resize((int(dim_ini[0]/red), int(dim_ini[1]/red))) ar = array(image_disque)# Définition de la matrice: tick_simu = 10 v = np.empty((duree_simu, np.shape(image_disque)[0], np.shape(image_disque)[1])) # conditions initiales: T_ini2 = 46.8 + 273.15 #K T_air2 = 25.3 + 273.15 #K T_metal = 35 + 273.15 #K #milieu du cercle: x_mi = np.shape(image_disque)[0]/2 y_mi = np.shape(image_disque)[1]/2 vide = np.zeros((np.shape(image_disque)[0], np.shape(image_disque)[1]), dtype=bool) for i in range(0, np.shape(image_disque)[0]): for j in range(0, np.shape(image_disque)[1]): pixel = image_disque.getpixel((j, i)) if pixel == (0, 247, 64, 255): v[0, i, j] = T_air2 vide[i, j] = True elif np.sqrt((i - x_mi)**2 + (j - y_mi)**2) >= 400 / red and np.sqrt((i - x_mi)**2 + (j - y_mi)**2) <= 448 / red and pixel != (0, 247, 64, 255): v[0, i, j] = T_ini2 else: v[0, i, j] = T_air2 dx = (16E-2) / np.shape(image_disque)[0] # largeur réelle divisée par nombre pixels A = (dt*lambd / ( rho * C_th * dx**2)) B = (2*dt*h) / (e*rho*C_th) C = (2*dt*epsilon * sigma ) / (e*rho*C_th) nX= np.shape(image_disque)[0] nY = np.shape(image_disque)[1] @njit def calculate(v): for k in range(1, tick_simu - 1): for i in range(1, nX-1): for j in range(0, nY-1): if vide[i, j] != True : v[k+1, i, j] = v[k, i, j] + A * (v[k, i+1, j] + v[k, i-1, j] + v[k, i, j+1] + v[k, i, j-1]) / 4 - (B * (v[k, i, j] - T_air2) + C * (v[k, i, j]**4 - T_air2**4)) return v v = calculate(v) plt.imshow(v[5], cmap='viridis', origin='lower', extent=(0, np.shape(image_disque)[1], 0, np.shape(image_disque)[0])) # Définissez les valeurs personnalisées pour l'axe des ordonnées et abscisses y_ticks = np.linspace(0, np.shape(image_disque)[0], 9) # Crée 9 valeurs de 0 à la hauteur de l'image x_ticks = np.linspace(0, np.shape(image_disque)[1], 9) # Crée 9 valeurs de 0 à la largeur de l'image plt.yticks(y_ticks, [-8, -6, -4, -2, 0, 2, 4, 6, 8]) # Personnalisez les étiquettes de l'axe des ordonnées plt.xticks(x_ticks, [-8, -6, -4, -2, 0, 2, 4, 6, '8 cm']) # Personnalisez les étiquettes de l'axe des abscisses plt.colorbar(label='Température (K)') plt.show()
j'ai corrigé "J'obtiens çaCode:calculate(v): for k in range(0, tick_simu - 1): for i in range(1, nX-1): for j in range(1, nY-1):hmmmm
Dernière modification par dembelekun ; 27/10/2023 à 10h36.
Dans la parenthèse derrière A*() il manque le morceau -4*v[k,i,j] que j'avais oublié.
Il faut tout initialiser : v[0, i, j] = T_air2 à remplacer par v[:, i, j] = T_air2
Mais le plus simple et le plus léger serait, à mon avis, d'avoir une matrice de température à 2 dimensions simplement
Sinon pour les CL sur les surfaces latérales, je vous envoie qqch.
Remarque : le temps caractéristique conductif est plutôt de 80 s (le 1300 s est pour le passage de la couronne vers le centre).
Dernière modification par gts2 ; 27/10/2023 à 13h11.
Comme ceci
Oui, mais il reste le pb
"imposer comme CL T=Tair sur les faces perpendiculaires au dessin, très douteux mais simple" est donc tellement douteux que cela ne marche pas.
J'ai un problème avec votre programme, sur mon ordinateur, lorsque je le lance j'ai l'image à t = 0s qui s'affiche assez rapidement mais lorsqu'il s'agit d'afficher une image après les 100s même au bout d'une heure rien ne s'affiche
Cela marche avec Spyder pas avec Pyzo !
Le problème est la fonction plt.show() de matplotlib est bloquante.
Il faut remplacer (ou supprimer ...) par plt.show(block=False) et peut être débuter par un plt.clf()
