r-热传递的球坐标边界条件的实现

Question

我想为半球应用热传递（热传导和对流）。它是球坐标系下的瞬态均匀传热。没有热量产生。半球的边界条件是在Tinitial=20度室温时开始的。外部环境温度为-22度。你可以想象半球是一种固体材料。此外，这是一个非线性模型，因为材料冻结后热导率会发生变化，这会改变温度分布。

我想找到这个固体在一定时间内的温度分布，直到中心温度达到-22度。

在这种情况下，温度取决于3个参数：T（r，theta，T）。半径、角度和时间。

1/α

我使用matlab应用了有限差分法，但是，边界条件存在问题。半球表面有对流，内部节点有传导，半球底部有恒定的温度，即空气温度（-22）。你可以在matlab文件中看到我为BC使用的脚本。

% 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;

程序运行后，新的温度值呈指数级增大，然后变为NaN。它应该向我展示固体的冷却和冷冻温度曲线，直到它达到塔尔温度。我想不出为什么会发生这样的变化。

我想听听你对BCs实施该计划的建议，或者我应该如何根据这种情况进行更改。提前感谢！！

Answer 1

您的代码太长，无法完全阅读和理解，但看起来您使用的是一个简单的前向Euler方案，对吗？如果是这样的话，试着减少时间步长dt，可能会减少很多，因为如果dt太大，这种方法可能会在数值上变得不稳定。这可能会降低计算速度（再次降低很多），但这就是你为这样一个简单的算法付出的代价。有一些替代方法不会出现不稳定性，但它们更难实现，因为你需要求解方程组。

很久以前，我用这个简单的方案做了一些热模拟。我发现稳定性标准是dt < (dx)^2 * c_p * rho / (6 * k)，这应该适用于在三维笛卡尔网格上的模拟，其中dx是空间步长，c_p是比热，rho是密度，k是材料的热导率。我不知道如何将其转化为球面坐标的情况。当时我学到的是选择小的时间步长，但更重要的是尽可能大的dx：当你将dx减少因子2时，你还需要将dt减少因子4以保持稳定。同时，对于3D问题，元素的数量将增加因子8。因此，1 / (dx)^5的总模拟时间尺度！！！