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