Mar-07-2017, 10:40 AM
Hi everyone. Thank you for reading.
We are working on solving this heat equation in a one dimension bar:
∂T/∂t=λ/ρc ∂²T/∂x²+(σ(T).i²)/(ρ.c.S²)
The second term (i²/...) is due to the Joule effect.
Notice that rho, lambda, sigma (electric conductivity) and c (the heat capacity) are function of temperature.
We have already tried to solve this equation with an implicite solving method but it actually doesn't work very well.
This scritp gives very weirds values (negative temperatures in Kelvins) and it's very (very) slow.
Does anyone have a solution to solve this type of equation? Thank you a lot for giving us your time!
Hugo.
We are working on solving this heat equation in a one dimension bar:
∂T/∂t=λ/ρc ∂²T/∂x²+(σ(T).i²)/(ρ.c.S²)
The second term (i²/...) is due to the Joule effect.
Notice that rho, lambda, sigma (electric conductivity) and c (the heat capacity) are function of temperature.
We have already tried to solve this equation with an implicite solving method but it actually doesn't work very well.
#liste des bibliotheques include import matplotlib.pyplot as plt import numpy as np from scipy.integrate import quad import mpmath as mm #effacement ecran avant simulation plt.clf() L=1 #1 metre T=1 #refroidissement sur T secondes D=1 #au pif.... M=50 #M points sur Ox, N intervalles de temps N=400 intensite=13000 largeur=0.02 #m hauteur=0.015 section=largeur*hauteur RRR=80 rho300K=17*10**-9 sigma0=rho300K/RRR alpha=3.93*10**-3 #coefficient temperature cuivre () mv=8960 tmax=500.0 #s idisplay=5#frequence d'affichage des courbes Ti=5 ### creation maillage x=np.linspace(0,L,N+1) ###calcul dx dx=x[1]-x[0] #calcul nb pas de temps dt=10.0 nbfourier=D*dt/dx**2 print ("nbfourier=",nbfourier) nbstep=int(tmax/dt)+1 print ('nbstep=',nbstep) print () #condition initiale TA=Ti*np.ones(N+1) TB=Ti*np.ones(N+1) #conductivité electrique def sigma(T): if T<=4.9: return 0 else : return sigma0*(1+alpha*T) #conductivité thermique def polynome1(T): return -98307.3620*T**4+159734.1213*T**3-62155.1495*T**2+21011.8166*T+4214.9608 def polynome2(T): return 13065.7087*T**4-118907.8785*T**3+403828.1132*T**2-606612.7583*T+340550.0003 def lamda(T): if T<4: return 0 #pour eviter les -infini if 4<=T<15: return polynome1(np.log10(T)) if 15<=T<=400: return polynome2(np.log10(T)) if T>400: return 380 ######################################################### Vmol=7.11*10**-6 largeur=0.02 #m hauteur=0.015 L=0.01 V=largeur*hauteur*L masse=mv*V n=masse/(mv*Vmol) kb=1.38064852*10**-23 Na=6.022140857*10**23 R=kb*Na γ=0.695 θD=493 def integrande(x): return (x**4*mm.exp(x)/((mm.exp(x)-1)**2)) def cp(T): return 9*n*R*(T/θD)**3*quad(integrande,0,θD/T)[0] ############################# def gamma(T): return dt*lamda(T)/(mv*cp(T)*dx**2) def beta(T): return dt*intensite**2*sigma(T)/(mv*cp(T)*section**2) #Initialisation du vecteur second membre BB=np.zeros(N+1) plt.clf() #affichage condition initiale plt.plot(x,TA) plt.show() #Boucle en temps for istep in range(1,nbstep): #Construction de la matrice AA=np.zeros((N+1,N+1)) ### x=0 AA[0,0]=1 ### 0 < x < L for i in range(1,N): AA[i,i-1]=-gamma(TA[i]) AA[i,i ]=1.+2.*gamma(TA[i]) AA[i,i+1]=-gamma(TA[i]) ### x=L AA[N,N]=1 RR=[beta(TA[i]) for i in range(N+1)] #Traitement des conditions aux limites #Calcul du vecteur second membre #------------------------------------- BB=TA.copy() CC=np.zeros(N+1) for i in range (N+1): CC[i]=BB[i]-RR[i] ####---Resolution du systeme lineaire TA=np.linalg.solve(AA,CC) #Affichage du champ de temperature dans barre if istep%idisplay==0: #afficher pas de temps print ("temps=",istep*dt) plt.xlim(0,L) plt.xlabel('$x$ [m]') plt.ylabel('$T$ [C]') plt.grid(True) plt.plot(x,TA,color=plt.get_cmap('jet')(float(istep)/nbstep)) plt.show() for i in range(0,N+1): TA[i]=TB[i] #--Dernier pas de temps plt.xlabel('$x$ [m]') plt.ylabel('$T$ [C]') plt.plot(x,TA,color=plt.get_cmap('jet')(0.999)) plt.show()(Sorry for the french terms)
This scritp gives very weirds values (negative temperatures in Kelvins) and it's very (very) slow.
Does anyone have a solution to solve this type of equation? Thank you a lot for giving us your time!
Hugo.