heat transfer for spherical coordinates boundary conditions implementation

Question

I want to apply heat transfer ( heat conduction and convection) for a hemisphere. It is a transient homogeneous heat transfer in spherical coordinates. There is no heat generation. Boundary conditions of hemisphere is in the beginning at Tinitial= 20 degree room temperature. External-enviromental temperature is -22 degree. You can imagine that hemisphere is a solid material. Also, it is a non-linear model, because thermal conductivity is changing after material is frozen, and this is going to change the temperature profile.

I want to find the temperature profile of this solid during a certain time until center temperature reach to -22 degree.

In this case, Temperature depends on 3 parameters : T(r,theta,t). radius, angle, and time.

1/α(∂T(r,θ,t))/∂t =1/r^2*∂/∂r(r^2(∂T(r,θ,t))/∂r)+ 1/(r^2*sinθ )∂/∂θ(sinθ(∂T(r,θ,t))/∂θ)

I applied finite difference method using matlab, However, boundary conditions have issues. There are convection on surface of the hemisphere, and conduction in the inner nodes, bottom of the hemisphere has constant temperature which is air temperature (-22). You can see the scripts which i am using for BCs in the matlab file.

% Temperature at surface of hemisphere solid boundary node

  for i=nodes
       for j=1:1:(nodes-1) 

Qcd_ot(i,j)=   ((k(i,j)+ k(i-1,j))/2)*A(i-1,j)*(( Told(i,j)-Told(i-1,j))/dr);         % heat conduction out of node

Qcv(i,j) =   h*(Tair-Told(i,j))*A(i,j); % heat transfer through convectioin on surface

Tnew(i,j)         =   ((Qcv(i,j)-Qcd_ot(i,j))/(mass(i,j)*cp(i,j))/2)*dt + Told(i,j);

      end       % end of for loop
     end

   % Temperature at inner nodes

   for i=2:1:(nodes-1)     
      for j=2:1:(nodes-1)  


 Qcd_in(i,j)=   ((k(i,j)+ k(i+1,j))/2)*A(i,j) *((2/R)*(( Told(i+1,j)-Told(i,j))/(2*dr)) + ((Told(i+1,j)-2*Told(i,j)+Told(i-1,j))/(dr^2)) + ((cot(y)/(R^2))*((Told(i,j+1)-Told(i,j-1))/(2*dy))) + (1/(R^2))*(Told(i,j+1)-2*Told(i,j)+ Told(i,j-1))/(dy^2));

 Qcd_out(i,j)=  ((k(i,j)+ k(i-1,j))/2)*A(i-1,j)*((2/R)*(( Told(i,j)-Told(i-1,j))/(2*dr)) +((Told(i+1,j)-2*Told(i,j)+Told(i-1,j))/(dr^2)) + ((cot(y)/(R^2))*((Told(i,j+1)-Told(i,j-1))/(2*dy))) + (1/(R^2))*(Told(i,j+1)-2*Told(i,j)+ Told(i,j-1))/(dy^2));


 Tnew(i,j)     =   ((Qcd_in(i,j)-Qcd_out(i,j))/(mass(i,j)*cp(i,j)))*dt + Told(i,j);



    end            %end  for loop
end                % end  for loop


   %Temperature for at center line nodes
  for i=2:1:(nodes-1)
      for j=1

     Qcd_line(i,j)=((k(i,j)+ k(i+1,j))/2)*A(i,j)*(Told(i+1,j)-Told(i,j))/dr;

     Qcd_lineout(i,j)=((k(i,j)+ k(i-1,j))/2)*A(i-1,j)*(Told(i,j)-Told(i-1,j))/dr;

     Tnew(i,j)= ((Qcd_line(i,j)-Qcd_lineout(i,j))/(mass(i,j)*cp(i,j)))*dt + Told(i,j);

      end 
 end


    % Temperature at bottom point (center) of the hemisphere solid
    for i=1
        for j=1:1:(nodes-1)

       Qcd_center(i,j)=(((k(i,j)+k(i+1,j))/2)*A(i,j)*(Told(i+1,j)-Tair)/dr);


       Tnew(i,j)= ((Qcd_center(i,j))/(mass(i,j)*cp(i,j)))*dt + Told(i,j);

      end
  end

   % Temperature at all bottom points of the hemisphere

     Tnew(:,nodes)=-22;


    Told=Tnew;

    t=t+dt;

Tnew temperatures values are getting bigger exponentially after program is run, and then becoming NaN. It supposed to show me cooling and freezing temperature profile of solid until it reaches to Tair temperature. I could not figure out the reasons why it is changing like that.

I would like to hear your suggestions for BCs implementation to this program, or how should i change them according to this conditions. Thanks in advance !!

Answer 1

Your code is too long to read and understand completely, but it looks like you are using a simple forward Euler scheme, is that correct? If so, try to reduce the time-step dt, maybe by a lot, since this method can become numerically unstable if dt is too big. This might slow down the speed of the computation (again by a lot), but that is the price you pay for such a simple algorithm. There are alternatives methods that do not suffer from instability, but they are much harder to implement, since you need to solve a system of equations.

I did some thermal simulations using this simple scheme a long time ago. I found that the stability criteria was dt < (dx)^2 * c_p * rho / (6 * k), which should be valid for a simulation on a 3D cartesian grid, where dx is the spatial step, c_p is the specific heat, rho the density and k the thermal conductivity of the material. I don't know how to convert this to your case with spherical coordinates. The thing I learned then was to choose small time-steps, but more importantly as large dx as possible: when you reduce dx by a factor 2, you also need to reduce dt by a factor 4 to keep things stable. At the same time, for a 3D problem, the number of elements will increase by a factor 8. So the total simulation time scales with 1 / (dx)^5!!!