#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ Created on Tue Mar 11 09:04:16 2025 @author: louisjeanmart """ #NYST chemical database import numpy as np from scipy.integrate import solve_ivp import matplotlib.pyplot as plt from scipy.optimize import newton from CoolProp.CoolProp import PropsSI from scipy.signal import savgol_filter import pandas as pd np.set_printoptions(threshold=np.inf) #------------------Database----------- # Charger le fichier Excel df = pd.read_excel('/Users/louisjeanmart/Desktop/Memoire/resultats_oxydation.xlsx') # Remplace par le bon chemin # Transformer chaque colonne en liste List_tga = [] for col in df.columns: List_tga.append(df[col].tolist()) Fe_540_Xtga = List_tga[0] Fe_540_Ytga = List_tga[4] Fe_540_Ytga2 = List_tga[5] def enthalpy_Fe(T): t = T/1000 # A = 23.97449 # B = 8.367750 # C = 0.000277 # D = -0.000086 # E = -0.000005 # F = 0.268027 # G = 62.06336 # H = 0 # return A * t + B * t**2 / 2 + C * t**3 / 3 + D * t**4 / 4 - E / t + F #- H if T >= 298 and T < 700: A = 18.42868 B = 24.64301 C = -8.913720 D = 9.664706 E = -0.012643 F = -6.573022 G = 42.51488 H = 0.000000 return A * t + B * t**2 / 2 + C * t**3 / 3 + D * t**4 / 4 - E / t + F #- H if T >= 700 and T <= 1042: A = -57767.65 B = 137919.7 C = -122773.2 D = 38682.42 E = 3993.080 F = 24078.67 G = -87364.01 H = 0 return A * t + B * t**2 / 2 + C * t**3 / 3 + D * t**4 / 4 - E / t + F #- H if T > 1042 and T < 1100: A = -325.8859 B = 28.92876 C = 0 D = 0 E = 411.9629 F = 745.8231 G = 241.8766 H = 0 return A * t + B * t**2 / 2 + C * t**3 / 3 + D * t**4 / 4 - E / t + F #- H if T >= 1100 and T < 1809: A = -776.7387 B = 919.4005 C = -383.7184 D = 57.08148 E = 242.1369 F = 697.6234 G = -558.3674 H = 0 return A * t + B * t**2 / 2 + C * t**3 / 3 + D * t**4 / 4 - E / t + F #- H def enthalpy_Fe0(T): t = T/1000 if T >= 298 and T < 1650: A = 45.75120 B = 18.78553 C = -5.952201 D = 0.852779 E = -0.081265 F = -286.7429 G = 110.3120 H = -272.0441 return A * t + B * t**2 / 2 + C * t**3 / 3 + D * t**4 / 4 - E / t + F #- H if T >= 1650 and T < 5000: A = 68.19920 B = -4.501232*1e-10 C = 1.195227*1e-10 D = -1.064302*1e-11 E = -3.092680*1e-10 F = -281.4326 G = 137.8377 H = -249.5321 return A * t + B * t**2 / 2 + C * t**3 / 3 + D * t**4 / 4 - E / t + F #- H def enthalpy_Fe3O4(T): t = T/1000 if T >= 298 and T < 900: A = 104.2096 B = 178.5108 C = 10.61510 D = 1.132534 E = -0.994202 F = -1163.336 G = 212.0585 H = -1120.894 return A * t + B * t**2 / 2 + C * t**3 / 3 + D * t**4 / 4 - E / t + F #- H if T >= 900 and T < 3000: A = 200.8320 B = 1.586435*10-7 C = -6.661682*10-8 D = 9.452452*10-9 E = 3.186020*10-8 F = -1174.135 G = 388.0790 H = -1120.894 return A * t + B * t**2 / 2 + C * t**3 / 3 + D * t**4 / 4 - E / t + F #- H def enthalpy_O2(T): t = T/1000 if T >= 100 and T < 700: A = 31.32234 B = -20.23531 C = 57.86644 D = -36.50624 E = -0.007374 F = -8.903471 G = 246.7945 H = 0 if T >= 700 and T < 2000: A = 30.03235 B = 8.772972 C = -3.988133 D = 0.788313 E = -0.741599 F = -11.32468 G = 236.1663 H = 0 if T >= 2000 and T < 6000: A = 20.91111 B = 10.72071 C = -2.020498 D = 0.146449 E = 9.245722 F = 5.337651 G = 237.6185 H = 0 return A * t + B * t**2 / 2 + C * t**3 / 3 + D * t**4 / 4 - E / t + F #- H Matrice_H = [] T_liste = np.arange(300,1800,1) #rafiner ca for i in range (len(T_liste)): h_liste = [T_liste[i],enthalpy_Fe(T_liste[i]),enthalpy_Fe0(T_liste[i]),enthalpy_Fe3O4(T_liste[i]),enthalpy_O2(T_liste[i])] #modifier la division par les masse mmolaire Matrice_H.append(h_liste) #------------------------------------- pi = np.pi R = 8.314 ####param initialisation Tg = 1088 #mettre a jour Temperature gas #1093 #Source papier k0_FeO = 2.67*1e-4 #grow rate parabolic FeO k0_Fe3O4 = 1.027*1e-6 #grow rate parabolic Fe3O4 Ta_FeO = 20319 #jsp FeO (je pense energie d'activation) Ta_Fe3O4 = 21310 #jsp Fe3O4 rho_FeO = 5745 #density Fe0 [kg/m3] rho_Fe3O4 = 5170 #density Fe3O4 [kg/m3] rho_Fe = 7874 #valeur intial rp_0 = 10*1e-6 #rp_0 = 31.5*1e-6 #pour poudre Mathieu X_0 = 10*1e-9 #X_FeO + X_Fe3O4 delta_0 = 1e-3 delta_Fe304_0 = 0.05 #expérience delta_Fe304 = X_Fe3O4/X_0 X_Fe3O4 = delta_Fe304_0*X_0 #initial couche oxide Fe3O4 X_FeO = X_0 - X_Fe3O4 #initial oxide couche FeO # vérifier r_Fe = rp_0 - X_0 A_FeO = 4 * pi * (r_Fe + X_FeO)**2 #Aire jusqu'a couche FeO Ap = 4 * pi * (r_Fe + X_FeO + X_Fe3O4)**2 #Aire particule #rapport massique stochio M_FeO = (55.845 + 15.999)/1000 M_O2 = (15.999*2)/1000 M_Fe = (55.845)/1000 M_FeO3 = 3*M_Fe + (4*15.999)/1000 nu_Fe_FeO = M_Fe/M_FeO nu_Fe_Fe3O4 = 3*M_Fe/M_FeO3 #### Second EDO bilan d'énergie #m_dot_R ????? nu_O2_FeO = M_O2/M_FeO #rapport massique O2/FeO stoich nu_O2_Fe3O4 = 2*M_O2/M_FeO3 #rapport massique O2/Fe2O3 stoich #radiative heat losses epsilon = 0.88 #emissivity sigma = 5.670374419 * 1e-8 #boltzman cst #Covective heat losses lambda_air = 6.763*1e-2 #considéré cst h_p = 2*lambda_air/(20*1e-6) #Re considéré = 0 air_elem = ['N2','O2'] air_mol = [0.79,0.21] #####Energie release q_Fe3O4 = 4.841*1e6 #source papers q_FeO = 3.787*1e6 h_O2 = enthalpy_O2(Tg)*1000 W_O2 = M_O2 #askip molar mass = molar weight pas en unité mais en valeur Matrice_h = [] def dHp_dt(dm_FeO_dt, dm_Fe3O4_dt, dm_O2_dt, Tp, Tg, Ap, q_FeO, q_Fe3O4, h_O2, W_O2, h_p, sigma, epsilon): term_chimique = q_FeO * dm_FeO_dt + q_Fe3O4 * dm_Fe3O4_dt - (h_O2 / W_O2) * dm_O2_dt term_convection = h_p * (Tp - Tg) term_radiation = sigma * epsilon * (Tp**4 - Tg**4) #term_radiation = 0 #negligeable return term_chimique - Ap * (term_convection + term_radiation) X_0 = X_FeO + X_Fe3O4 r_Fe = rp_0 - X_0 A_FeO = 4 * np.pi * (r_Fe)**2 Ap_Fe3O4 = 4 * np.pi * (r_Fe + X_FeO)**2 mFeO_0 = (rho_FeO * A_FeO * X_FeO) m_Fe3O4_0 = (rho_Fe3O4 * Ap_Fe3O4 * X_Fe3O4) m_Fe_0 = 4/3 * np.pi * r_Fe**3 * rho_Fe m0 = m_Fe_0 + m_Fe3O4_0 + mFeO_0 m0X = m0/M_Fe*(M_FeO3/3 + M_FeO) Tp = 300 H0 = (m_Fe_0/ M_Fe*Matrice_H[Tp-300][1]+mFeO_0/M_FeO*Matrice_H[Tp-300][2]+m_Fe3O4_0/M_FeO3*Matrice_H[Tp-300][3])*1000 test = [] X=np.zeros(2) X1 = [] X2 = [] Rayon= [] def reaction_odes(t, m, params,isoT,Tg): # Initialisation des valeurs initiales (nu_Fe_FeO, nu_Fe_Fe3O4, rho_FeO, k0_FeO, Ta_FeO, rho_Fe3O4, k0_Fe3O4, Ta_Fe3O4) = params if isoT == True: T_guess = 300 else: if t==0: T_guess = 300 else: T_guess = test[-1] Tp = T_guess r_Fe = (m[2]/(rho_Fe*4/3*pi))** (1/3) #me parait bon mais a voir Rayon.append(r_Fe) #print(m[2]) #r = rp_0 + X_more X[0] = m[0]/(4*pi*(r_Fe**2)*rho_FeO) X[1] = m[1]/(4*pi*(r_Fe+X[0])**2*rho_Fe3O4) A_FeO = 4 * np.pi * (r_Fe + X[0])**2 Ap = 4 * np.pi * (r_Fe + X[0] + X[1])**2 X1.append(X[0]) X2.append(X[1]) dm_FeO_dt = (rho_FeO * A_FeO * k0_FeO / X[0]) * np.exp(-Ta_FeO / Tp) dm_Fe3O4_dt = (rho_Fe3O4 * Ap * k0_Fe3O4 / X[1]) * np.exp(-Ta_Fe3O4 / Tp) dm_Fe_dt = -nu_Fe_FeO * dm_FeO_dt - nu_Fe_Fe3O4 * dm_Fe3O4_dt dm_O2_dt = -nu_O2_FeO * dm_FeO_dt - nu_O2_Fe3O4 * dm_Fe3O4_dt Tsf = (2*Tp + Tg)/3 lambda_air = 0 for i in range(len(air_elem)): lambda_air += (PropsSI('L', 'T', Tsf, 'P', 101325, air_elem[i])*air_mol[i]) lambda_air += 1/(air_mol[i]/PropsSI('L', 'T', Tsf, 'P', 101325, air_elem[i])) lambda_air *=0.5 h_p = 2*lambda_air/(20*1e-6) #Re considéré = 0 Hp = dHp_dt(dm_FeO_dt, dm_Fe3O4_dt, dm_O2_dt, Tp, Tg, Ap, q_FeO, q_Fe3O4, h_O2, W_O2, h_p, sigma, epsilon) l = [] z = 0 mem = 0 for i in range(len(Matrice_H)): Hp_new = ((m[2])/M_Fe*Matrice_H[i][1] + (m[0]/M_FeO)*Matrice_H[i][2] + (m[1]/M_FeO3)*Matrice_H[i][3] )*1000 l.append(Hp_new) z = 10000 # Initialiser avec une valeur infinie T_test = Tp for i in range(len(l)): diff = abs(l[i] - m[4]) # Calculer la différence entre H_new et H0 if diff < z: z = diff mem = i T_test = Matrice_H[i][0] T_guess = Tp Tp = T_test Matrice_h.append(l[-1]) test.append(Tp) return [dm_FeO_dt, dm_Fe3O4_dt,dm_Fe_dt,dm_O2_dt,Hp] #return [X[0], X[1]] #detailler les courbes et tjrs voir le problème de temperature def euler_explicite(f, y0, t0, tf, h, params,Tg): n = int((tf - t0) / h) + 1 t = np.linspace(t0, tf, n) y = np.zeros((n, len(y0))) y[0] = y0 for i in range(n - 1): dydt = f(t[i], y[i], params,False,Tg) y[i + 1] = y[i] + h * np.array(dydt) return t, y # Paramètres d'exemple params = (nu_Fe_FeO, nu_Fe_Fe3O4, rho_FeO, k0_FeO, Ta_FeO, rho_Fe3O4, k0_Fe3O4, Ta_Fe3O4) m_0 = [mFeO_0, m_Fe3O4_0,m_Fe_0,0,H0] t_span = (0, 0.03) h = 0.0001 t_list = [] y_list =[] temp = [] t, y = euler_explicite(reaction_odes, m_0, t_span[0], t_span[1], h, params,Tg) #Qd pas thermal runaway # y_smooth = savgol_filter(test, window_length=1001, polyorder=5) #permet de smooth la fonction et eviter les problèmes due a la matrice # y_smooth = savgol_filter(test, window_length=101, polyorder=8) # T_liste = [1082,1090,1092,1093] # for i in range(len(T_liste)): # test = [] # t, z = euler_explicite(reaction_odes, m_0, t_span[0], t_span[1], h, params,T_liste[i]) # t_list.append(t) # y_list.append(z) # temp.append(test) # plt.plot(savgol_filter(t_list[3][:-1], window_length=101, polyorder=8), temp[3], label='r') # plt.plot(savgol_filter(t_list[0][:-1], window_length=1001, polyorder=5), temp[0], label='r') # plt.plot(savgol_filter(t_list[1][:-1], window_length=1001, polyorder=5), temp[1], label='r') # plt.plot(savgol_filter(t_list[2][:-1], window_length=1001, polyorder=5), temp[2], label='r') # plt.plot([t_list[3][548], t_list[3][548]], [1100, temp[3][549]], linestyle='--', color='gray') # ligne verticale # plt.plot([t_list[3][0], t_list[3][548]], [temp[3][549], temp[3][549]], linestyle='--', color='gray') # ligne horizontale # plt.xlim(0.0035,0.0099) # plt.ylim(1100,1600) # plt.xticks([]) # plt.yticks([]) # ax = plt.gca() # ax.spines['top'].set_visible(False) # ax.spines['right'].set_visible(False) # #plt.xlabel('Time (s)') # plt.savefig("Termal_runaway.pdf") # #plt.ylabel('Temperature ') # plt.show() plt.plot(t[:-1], test, label='T') #plt.plot(t[:-1], y_smooth, label='T') plt.show()