Heat equation

[HugoL]

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.

#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.

[micseydel]

Since your script is very long, and the question is more about physics than Python, you might have more luck collaborating with classmates. Someone might come along here, but in the event that doesn't happen, I didn't want to leave you hanging.