# Heat equation problem

**URL:** <https://community.freefem.org/t/heat-equation-problem/386>\
**Category:** General Discussion\
**Created:** [May 4, 2020, 12:56pm UTC](https://community.freefem.org/t/heat-equation-problem/386 "2020-05-04T12:56:06Z")\
**Posts on this page:** 1\
**Page:** 1

<div class="post-metadata">

**Author:** ![leopold](https://avatars.discourse-cdn.com/v4/letter/l/6bbea6/32.png) [@leopold](https://community.freefem.org/u/leopold)\
**Post date:** [May 4, 2020, 12:56pm UTC](https://community.freefem.org/t/heat-equation-problem/386/1 "2020-05-04T12:56:06Z")

</div>

Hello,

I am trying to solve the heat equation in a molten salt and in a metallic tank in 1D. The thickness of the system is equal to 7cm (2cm for the tank wall and 5cm for the salt). The salt has a residual power and initially its temperature is equal to 1050°C. The initial temperature of the tank is 20°C. There is radiation at the top of the salt and convection at the bottom of the tank (optional). I supposed that there is conduction in the whole system. The temperature in the air at the top and at the bottom is equal to 20°C.

I wrote the code in Python with finite difference method to compare the results and they are completely different. Here is the FreeFem code :

int np=70; //nombre de noeuds  
mesh Th = square(np,1,[0.07 \* x,y/1000]); //maillage

//plot(Th,wait=1,aspectratio=1);

real Trecup=20+273, TselInit=1050+273;  
func rho = 8860 \* (x\<=0.02)+1527 \* (x\>0.02);  
func cp = 586 \* (x\<=0.02)+1155 \* (x\>0.02);  
func kth = 23.6 \* (x\<=0.02)+0.55 \* (x\>0.02);  
func u0 = Trecup \* (x\<=0.02)+TselInit \* (x\>0.02);  
real rad=1e-8, emissivite=0.2, dt=1, duree=1000, hg=0, uSup=20+273, uInf=20+273, pInit=3e+9;  
func real fsource(real t){ //calcul puissance résiduelle (W/m3)  
if(t\<=1){  
return (0.0377 \* pInit/18);}  
if((t\>1) && (t\<=100)){  
return ((-0.384 \* log(t)+3.77)/100) \* pInit/18;}  
if((t\>100) && (t\<=3600)){  
return (-0.279 \* log(t)+3.285)/100 \* pInit/18;}  
if((t\>3600) && (t\<=86400)){  
return (-0.189 \* log(t)+2.546)/100 \* pInit/18;}  
if((t\>86400) && (t\<=2.628E6)){  
return (-0.0673 \* log(t)+1.165)/100 \* pInit/18;}  
if((t\>2.628E6) && (t\<=3.154E7)){  
return (-0.0624 \* log(t)+1.093)/100 \* pInit/18;}  
if((t\>3.154E7) && (t\<=1E9)){  
return (-0.00304 \* log(t)+6.743E-2)/100 \* pInit/18;}  
if(t\>1E9){  
return 0;}  
}  
real source=fsource(0);

fespace Vh(Th,P1);  
Vh vold,w,v=u0,b;

problem thermradia(v,w)= int2d(Th)(rho \* cp \* v \* w/dt + kth \* (dx(v) \* dx(w) + dy(v) \* dy(w)))+int1d(Th,2)(b \* w)- int2d(Th)(rho \* cp \* vold \* w/dt)- int2d(Th)(source \* (x\>0.02) \* w)+ int1d(Th,4)(hg \* v \* w)- int1d(Th,4)(hg \* uInf \* w);

for(real t=0;t\<duree;t+=dt){  
cout\<\<t\<\<endl;  
plot(v,value=true,fill=true);  
vold=v;  
source = fsource(t);  
for(int m=0;m\<5;m++){  
b= rad \* emissivite \* (v - uSup) \* (v + uSup) \* ((v+uSup) \* (v+uSup) + uSup \* uSup);  
thermradia;  
}  
}  
plot(v,value=true,fill=true,dim=3,ps=“1DFreeFem.eps”);

Is there something wrong in my code ?  
Thanks in advance,

Leopold