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

Question

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

我想在中心温度达到-22度之前找到这种固体的温度曲线。

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

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

我使用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;

Tnew 温度值在程序运行后呈指数增长，然后变为 NaN。它应该向我显示固体的冷却和冷冻温度曲线，直到它达到 Tair 温度。我无法弄清楚它为什么会这样变化的原因。

我想听听您对本程序实施 BC 的建议，或者我应该如何根据这种情况更改它们。提前致谢！！

Answer 1

您的代码太长，无法完全阅读和理解，但看起来您使用的是简单的forward Euler 方案，对吗？如果是这样，请尝试减少时间步长dt，可能会减少很多，因为如果dt 太大，此方法可能会在数字上变为unstable。这可能会减慢计算的速度（再次减慢很多），但这就是您为这样一个简单的算法付出的代价。有alternatives methods 不会受到不稳定的影响，但它们实现起来要困难得多，因为您需要求解方程组。

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