Thermodynamique : équation de diffusion

[dembelekun]

Bonjour,


Dans le cadre d'un devoir, je m'intéresse au refroidissement d'un disque de frein en prenant en compte le transfert thermique par convection et rayonnement, seulement je n'ai pas encore vu les laplacien et donc je galère un peu pour établir l'équation.
Le but à la fin, est d'avoir une EDP et la résoudre sur python pour avoir une jolie représentation graphique. Dans ce genre : https://levelup.gitconnected.com/sol...n-3334004aa01a


Je ne sais pas si j'ai intérêt à utiliser le laplacien pour les repères cartésien ou celui pour les repères cylindriques.

J'ai :

transfert par rayonnement :
transfert par convection :

Puis l'équation de diffusion, :


La question est donc, comment relier P reçue avec ces deux transferts thermiques ?
-----

[gts2]

Je ne sais pas si j'ai intérêt à utiliser le laplacien pour les repères cartésien ou celui pour les repères cylindriques.

Le disque étant cylindrique, la réponse en analytique serait évidemment cylindrique (disparition du \theta), en numérique il y a peut-être du pour et du contre, en particulier gestion du centre r=0.

La question est donc, comment relier P reçue avec ces deux transferts thermiques ?

Votre P reçue, j'ai plutôt tendance à l'appeler produite ou créée, de toute manière dans un frein, elle est nulle.

[gts2]

Précisions :

Il manque le z dans le laplacien.

Les deux transferts thermiques sont des conditions aux limites en z=0 z=h et r=R.

[dembelekun]

Je ne peux pas uniquement m'intéresser à la température de surface ? et donc ne pas prendre en compte l'épaisseur

[A voir en vidéo sur Futura]

[dembelekun]

très bien ,qu'est ce qui va changer mis à part le laplacien ?

[gts2]

Je ne peux pas uniquement m'intéresser à la température de surface ? et donc ne pas prendre en compte l'épaisseur.

Sous l'hypothèse système mince (température intérieure uniforme selon z à justifier), oui, et le plus simple, plutôt que d'utiliser des équations toutes faites, est d'écrire le premier principe sur une couronne d'épaisseur e entre r et r+dr, cela donne une équation aux dérivées partielles en T(r,t).

[dembelekun]

J'ai retrouvé ça sur mon ordinateur, j'avais ça il y a longtemps. Seulement, je ne suis pas vraiment sûr

[gts2]

C'est bien l'idée, mis il y a manifestement un certain nombre de fautes de frappe :
: oubli de dT, un d \theta qui est un dr ...
: 1/r à la place de e ...
: 1/r oublié ...

[dembelekun]

Merci beaucoup, vous vous y connaissez en Cranck-Nilcoson cylindrique ? J'ai fait quelques recherches ça me semble barbare

[gts2]

Quel est le but du jeu ?

Une méthode d'Euler simple fonctionne relativement bien (il faut juste faire attention à la stabilité, donc prendre un pas temporel suffisamment petit)

[dembelekun]

Pour résoudre l'équation de diffusion, comme dit au début du poste, j'aimerais bien qu'à la fin, graphiquement, j'ai la même chose que l'article que j'ai envoyé

[gts2]

Justement, l'article cité utilise une méthode d'Euler basique.

[dembelekun]

Ce qui nous donne ça ?

[gts2]

Vous avez oublié le r (équation inhomogène)

Je trouve plutôt

[dembelekun]

Ok pour le r et l'erreur de signe mais vous n'avez pas de 2dr^2 ?

[gts2]

Le deuxième terme a bien un /2 : différence centrée.

[dembelekun]

très bien, j'ai donc cette suite Capture d’écran 2023-08-11 173346.jpg

(le code latex si vous souhaitez modifier quelque chose: u_{r}^{k+1} = \Delta t \frac{\lambda}{r \rho c_{th}} (\frac{u_{r+1}^{k}-u_{r-1}^{k}}{2\Delta r} + r\frac{u_{r+1}^{k}+u_{r-1}^{k} - 2u_{r}^{k}}{\Delta r^2}) - \frac{\Delta t }{e\rho C_{th}}(h(u_{r}^{k} - T_{\infty}) + \epsilon \sigma (u_{r}^{k}^{4}-T^{4}_{\infty})) )


J'essaye depuis un bon moment de coder ça sur python mais je ne comprends pas, dès la première seconde les températures tombent à 0 Kelvin. Je n'arrive pas à trouver d'où ça vient Figure 2023-08-11 173100.png

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from matplotlib.animation import FuncAnimation

plate_length = 50
max_iter_time = 10000

dt = 0.01
dr = 1

#rayon
rayon = []
for i in range(0,50):
    rayon.append(i)
    
#paramètre
T_ini = 46.8 + 273.15 #K
T_air = 25.3 + 273.15 #K 
C_th =  502 #J.kg-1.K-1
e = 0.003 #m 
rho = 7850 # kg.m-3
h = 20 #W.m-2.K-1
lambd = 20 #W.m-1.K-1
epsilon = 0.15
sigma = 5.67E-8 #W.m-2.K-4


#Initialize solution: the grid of u(k, r)
u = np.empty((max_iter_time, plate_length))

#Conditions initiales:
u_initial = T_ini 
u.fill(u_initial)

    

#conditions limites
u[:,:20] = T_air
u[0,30:] = T_air
u[:,20:30] = T_ini

#Calculs EDP
def calculate(u):
    for k in range(0, max_iter_time - 1):
        for r in range(1, plate_length - 1):
            u[k + 1, r] = (dt*lambd / ( r*rho * C_th)) * (( (u[k, r + 1] - u[k, r - 1]) / (2 * dr)) + r * (u[k, r + 1] + u[k, r - 1] 
            - 2 * u[k, r]) /  (2 * dr**2)) -(dt/(e * rho * C_th)) * (h * (u[k, r] - T_air) + epsilon * sigma * (u[k, r] ** 4 - T_air ** 4))

    return (u)
    
u = calculate(u)

# Plot the temperature distribution at a specific time instant
time_instant = 100  # Choose the time instant you want to plot

plt.plot(rayon, u[time_instant, :])
plt.xlabel('Radius')
plt.ylabel('Temperature (K)')
plt.title(f'Temperature Distribution at t = {time_instant * dt:.2f} seconds')
plt.grid(True)
plt.show()

[dembelekun]

Je viens de trouver le problème, j'ai oublié le à la fin

[gts2]

Il y a d'autres problèmes :

Dans la boucle r est l'index repérant le rayon pas le rayon, donc le r qui intervient en tant que variable (pas indice) est qqch du genre rR=r*rayon/plate_length, rayon = rayon réel en mètre et rR r réel.
Idem pour le dr.

J'ai du mal avec le slicing en Python, donc prendre mes remarques avec des bémols :
#conditions limites
u[:,:20] = T_air n'a aucun effet puisque vous recalculez u[k,x] pour x allant de 1 à r max
u[0,30:] = T_air
u[:,20:30] = T_ini pour moi la valeur initiale est pour u[0,...] donc u[0,:] = T_ini ?
On a l'impression que le frein est entre r=20 et r=30 ?

La condition aux limites Tair est trop forte, cela impose que le frein au niveau de votre limite prenne instantanément la température de l'air.
Il vaudrait mieux considérer au niveau des faces internes et externes que le flux est donné par votre expression.

[dembelekun]

Merci pour vos remarques,
J'ai modifié pour le rayon mais est-ce nécessaire de toucher à dr ?

Effectivement, les cases qui définissent la zone du disque de frein ont été choisi arbitrairement, je modifierai ça plus tard, c'était histoire de voir si ça fonctionnait bien.
Quant à la phrase "Il vaudrait mieux considérer au niveau des faces internes et externes que le flux est donné par votre expression." je ne pense pas avoir tout compris. Encore une fois, merci pour votre aide


J'ai un petit problème avec ma température pour les rayons proches de 0, les températures augmentent grandement

 
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from matplotlib.animation import FuncAnimation

plate_length = 50
max_iter_time = 1000

dt = 1
dr = 0.01

#rayon
rayon = []
for i in range(0,50):
    rayon.append(i)
    
#paramètre:
T_ini = 46.8 + 273.15 #K
T_air = 25.3 + 273.15 #K 
C_th =  502 #J.kg-1.K-1
e = 0.003 #m 
rho = 7850 # kg.m-3
h = 20 #W.m-2.K-1
lambd = 20 #W.m-1.K-1
epsilon = 0.15
sigma = 5.67E-8 #W.m-2.K-4
r_disc = 8E-2 #cm 

m = r_disc / plate_length 


#Initialize solution: the grid of u(k, r)
u = np.empty((max_iter_time, plate_length))

#Conditions initiales:
u_initial = T_ini 
u.fill(u_initial)

    

#conditions limites:
u[:,:10] = T_air
u[:,40:] = T_air

u[0,10:40] = T_ini

#Calculs EDP:
def calculate(u):
    for k in range(0, max_iter_time - 1):
        for r in range(1, plate_length - 1):
            u[k + 1, r] = (dt*lambd / ( m*r*rho * C_th)) * (( (u[k, r + 1] - u[k, r - 1]) / (2 * dr)) + r * m * (u[k, r + 1] + u[k, r - 1] 
            - 2 * u[k, r]) /  (2 * dr**2)) -(dt/(e * rho * C_th)) * (h * (u[k, r] - T_air) + epsilon * sigma * (u[k, r] ** 4 - T_air ** 4)) +u[k,r]

    return (u)
    
u = calculate(u)


# #Tracé animation:
# def plotheatmap(u_k, k):
#     # Clear the current plot figure
#     plt.clf()
    
#     plt.plot(rayon, u[k, :])
#     plt.title(f"Temperature at t = {k*dt//60:.0f} min" f"{k*dt%60:.0f} sec")
#     plt.xlabel("x")
#     plt.ylabel("Température(K)")
#     plt.ylim(298,330)
    
    
    
    
#     return plt

# # Do the calculation here
# u = calculate(u)

# def animate(k):
#     plotheatmap(u[k], k)

# anim = animation.FuncAnimation(plt.figure(), animate, interval=1, frames=max_iter_time, repeat=False)
# anim.save("heat_equation_solution.gif")

# print("Done!")




# Calculer les valeurs de temps et de position pour l'axe des abscisses
num_levels = 15  # Nombre de niveaux de couleurs
temp_levels = np.linspace(np.min(u), np.max(u), num_levels)

# Créer la figure et le graphique
fig, ax = plt.subplots()
contour = ax.contourf(positions_grid, times_grid, u.T, levels=temp_levels, cmap='viridis') # ou cmap = 'hot'
cbar = plt.colorbar(contour, ax=ax, label='Température (K)')
ax.set_xlabel('Position (m)')
ax.set_ylabel('Temps (s)')
ax.set_title('Distribution de la température dans le disque de frein')

# Afficher le graphique
plt.show()

[gts2]

Je regarde de plus près mais les données sont incohérentes :
on a un dr de 0,01 et vous faites variez l'indice de 0 à 50 donc r de 0 à 50 cm et vous annoncez un rayon du disque de 8 cm.
ensuite vous imposez Tair de 0 à 10 donc de 0 à 10 cm (trou central ?) et de 40 à 50 donc de 40 cm à 50 cm, cela voudrait dire que le disque a un rayon externe de 40 cm ? et un trou interne de 10 cm ? pas plus cohérent avec vos 8 cm
vous imposez des conditions aux limites de 0 et 10 et de 40 à 50 mais vous faites quand même le calcul de T de 1 à 49.

[dembelekun]

Alors voilà mon idée, les cases de [1;49] seront dédiées à la zone du disque de frein (environ 8 cm). Les colonnes 0 et 50 seront mon thermostat (l'air.) Le patin ne recouvrant pas l'intégralité de mon disque de frein, seulement une partie du disque sera à la température initiale.

Certes le disque n'est pas plein, mais les tiges qui relie la périphérie (?) au moyeu (?) aide à la dissipation et la chaleur s'y propage

[gts2]

Alors voilà mon idée, les cases de [1;49] seront dédiées à la zone du disque de frein (environ 8 cm). Les colonnes 0 et 50 seront mon thermostat (l'air.) Le patin ne recouvrant pas l'intégralité de mon disque de frein, seulement une partie du disque sera à la température initiale.

Dans ce cas, à quoi correspond :
#conditions limites:
u[:,:10] = T_air
u[:,40:] = T_air
u[0,10:40] = T_ini

Si le frein est de 8cm que vous découpez en 50, dr=8cm/50=1,6 mm et non dr=1cm.

Certes le disque n'est pas plein, mais les tiges qui relie la périphérie (?) au moyeu (?) aide à la dissipation et la chaleur s'y propage

Si vous considérez le disque comme plein, en considérant que le moyeu se comporte comme le disque, il va y avoir en effet des problèmes en r=0 (le 1/r du laplacien), mais dans ce cas il n'y a plus de conditions aux limites T=Tair en r=0.

SI vous faisiez un dessin cela serait plus clair.

Vous ne pouvez prendre comme condition aux limites T=Tair, cela est incohérent avec votre modèle qui repose sur de la convection qui suppose justement Tsurface différent de Tambiant.

[gts2]

Pour la condition aux limites du rayon extérieur, on écrit la conservation du flux que l'on discrétise , on en déduit alors T(R).

[gts2]

Bonjour,

Un disque avec le "moyeu" ressemble à cela (wikipedia) :

Donc la partie centrale est nettement différente, vous ne pouvez appliquer la même loi thermique.
Une possibilité : considérer que c'est une "grosse" pièce donc un thermostat à température Tair (ce que vous faites avec u[:,:10] = T_air), mais dans ce cas c'est un thermostat et il ne faut pas lui appliquer la loi thermique :
for r in range(1, plate_length - 1): devrait être for r in range(10, plate_length - 1):

[dembelekun]

J'ai un disque de frein de vélo, j'aurais du le préciser. J'ai considéré 3 zone, le centre allant de 0 à 2cm, je le considère comme vide. De 2 à 6cm, les branches qui ne sont pas en contact avec le patin donc à t = 0 sont à la température de l'air et enfin de 6 à 7,5cm la zone chauffée. J'ai laissé un thermostat à l'extérieur du disque disque-de-frein-vtt-160180203mm-cooma-wavy-mtb-rotor-wavy-160180203-neuf.jpg

J'ai modifié ma boucle for sur r pour ne pas avoir à recalculé les températures du thermostat. Cependant, j'ai des problèmes sur les extrémités, votre message sur les conditions limites est-il la solution ?
Figure 2023-08-13 103433.png

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from matplotlib.animation import FuncAnimation
from matplotlib.colors import LinearSegmentedColormap
plate_length = 50
max_iter_time = 1000

dt = 1
dr = 0.01

#rayon
rayon = []
for i in range(0,50):
    rayon.append(i)
    
#paramètre:
T_ini = 46.8 + 273.15 #K
T_air = 25.3 + 273.15 #K 
C_th =  502 #J.kg-1.K-1
e = 0.003 #m 
rho = 7850 # kg.m-3
h = 20 #W.m-2.K-1
lambd = 20 #W.m-1.K-1
epsilon = 0.15
sigma = 5.67E-8 #W.m-2.K-4
r_disc = 8E-2 #cm 

m = r_disc / plate_length 


#Initialize solution: the grid of u(k, r)
u = np.empty((max_iter_time, plate_length))

#Conditions initiales:
u_initial = T_ini 
u.fill(u_initial)

    

#conditions limites:
u[:,:12] = T_air # zone interieure du disque thermostat
u[:,48:] = T_air # zone exterieure du disque thermostat

u[0,35:47] = T_ini #bande du disque chauffé
u[0,1:34] = T_air # centre du disque

#Calculs EDP:
def calculate(u):
    for k in range(0, max_iter_time - 1):
        for r in range(13,  47 ):
            u[k + 1, r] = (dt*lambd / ( m*r*rho * C_th)) * (( (u[k, r + 1] - u[k, r - 1]) / (2 * dr)) + r * m * (u[k, r + 1] + u[k, r - 1] 
            - 2 * u[k, r]) /  ( dr**2)) -(dt/(e * rho * C_th)) * (h * (u[k, r] - T_air) + epsilon * sigma * (u[k, r] ** 4 - T_air ** 4)) +u[k,r]

    return (u)
    
u = calculate(u)


# Calculer les valeurs de temps et de position pour l'axe des abscisses
num_levels = 10  # Nombre de niveaux de couleurs
temp_levels = np.linspace(np.min(u), np.max(u), num_levels)
times_grid, positions_grid = np.meshgrid(np.arange(max_iter_time), np.arange(plate_length))

# Créer la figure et le graphique
fig, ax = plt.subplots()
contour = ax.contourf(positions_grid, times_grid, u.T, levels=temp_levels, cmap='RdYlBu_r')
cbar = plt.colorbar(contour, ax=ax, label='Température (K)')

# Transformation linéaire des valeurs sur l'axe des x
original_x_values = np.linspace(0, 50, 5)  # Valeurs d'origine de 0 à 50
transformed_x_values = [0,2,4,6,8]  # Valeurs transformées de 0 à 8

# Associer les positions des marques d'échelle avec les valeurs transformées
ax.set_xticks(original_x_values)
# Utiliser les valeurs d'origine pour les étiquettes, mais afficher les valeurs transformées
ax.set_xticklabels(transformed_x_values)

ax.set_xlabel('Position (cm)')
ax.set_ylabel('Temps (s)')
ax.set_title('Distribution de la température dans le disque de frein')

# Afficher le graphique
plt.show()

[dembelekun]

ok, en mettant for r in range(12, 48 ): ça à l'air OK. J'aimerais avoir votre avis

[gts2]

Je pense plutôt à un problème de slicing :

Je raisonne avec des paramètres (dès qu'une même valeur apparait plus d'une fois, il faut lui donner un nom)
r1=12
u[:,:r1] = T_air # zone interieure du disque thermostat, remplit de 0 à r1-1
r3=48
u[:,r3:] = T_air # zone exterieure du disque thermostat, remplit de r3 à ...
r2=35
u[0,r2:r3] = T_ini #bande du disque chauffé, remplit de r2 à r3-1 (autrement dit mettre 48 et pas 47)
u[0,0:r2] = T_air # centre du disque, remplit de 0 à r2-1 (autrement dit mettre 35 et pas 34)

boucle
for r in range(r1, r3): boucle de r1 à r3-1 (autrement dit 12,48

Autrement dit, il faut vous habituer au fait que les informaticiens comptent de 0 à N-1 alors que le reste du monde compte de 1 à N (je n'ai jamais vraiment pu m'y habituer)

Autre remarque d'optimisation cette fois, vous calculez 50*1000 fois dt*lambd / ( rho * C_th) ; il faut le calculer une fois pour toutes avec a1=dt*lambd / ( rho * C_th* dr**2) (autant inclure le dr une fois pour toutes).

Dernière remarque dr et m c'est la même chose : variation de r correspondant à une augmentation d'indice de 1.

[dembelekun]

Très bien, je modifierai ça. Merci infiniment pour votre aide

[dembelekun]

Bonsoir, je vous recontacte puisque j'ai du nouveau. Mes professeurs m'ont dit que ça serait pas mal de prendre en compte la géométrie du disque dans son intégralité. J'ai donc fait un petit programme python me permettant d'avoir la forme du disque

Je souhaiterai savoir si je peux reprendre le programme précédent ( resté en cylindrique ) ou alors je vais devoir passer en cartésien ? Si je reste en cylindrique comment modifier mon programme ?

Je vous remercie pour votre aide