使用matlab进行球坐标的传热答案

【问题标题】：heat transfer for spherical coordinates using matlab使用matlab进行球坐标的传热
【发布时间】：2013-04-09 15:25:19
【问题描述】：

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

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

在这种情况下，温度取决于 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 应用了有限差分法，但是程序不会为半球的内部节点计算任何东西，而只是给我初始温度值（这里告诉了）。您可以看到我用于内部节点的一些脚本。

% initial conditions
Tair    =   -30.0;               % Temperature of air
Tin     =    21; 

% setting initial values for grid
for i=1:(nodes)
    for j=1:(nodes)
Told(i,j)     =   Tin;
Tnew(i,j)     =   Tin;
frozen(i)   =   0;                      
latent(i)   =   Qs*mass(i)*Water/dt;    
k(i)        =   0.5;                    
cp(i)       =   cw;                     
W(i)        =   Water;                  
l(i)        =   0;                      
S(i)        =   1-Water;                
     end
end


%Simulation conditions
J       =    9;              % No. of space steps   
nodes   =   J+1;             % Number of nodes along radius or theta direction
dt      =0.1;


t       =   0;               % time index on start
tmax    =   7000;            % Time simmulation is supposed to run


R       =  d/2;

dr      =  (d/2)/J;         %  space steps in r direction

 y   =        pi/2;        % (theta) for hemisphere
dy      =  (pi/2)/J;       % space steps in Theta direction

    % Top surface condition for hemisphere
      i=nodes; 
        for j=1:1:(nodes-1) 

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

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

        Tnew(i,j) = ((Qcv(i,j)-Qcd_ot(i,j))/(mass(i)*cp(i))/2)*dt + Told(i,j); 
          end             %end of for loop 

  % Temperature profile for inner nodes

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

Qcd_in(i,j)=   ((k(i)+ k(i+1))/2)*A(i) *((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))/(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)+ k(i-1))/2)*A(i-1)*((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)-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)*cp(i)))*dt + Told(i,j);
    end           
end               

   %bottom of the hemisphere solid
  Tnew(:,nodes)=-30;



 Told=Tnew;
 t=t+dt;

编辑 *谢谢，现在脚本正在运行和计算。我可以看到模型系统的温度曲线。

但是，我想为此半球温度曲线绘制 2D 或 3D 图。另外，如果可能的话，我想在特定时间运行温度变化动画。我用于模拟和编写文件的代码是

     t=0;

    tmax=7000;
    ...................
    .....................
    ss=0;   % index for printouts

     %start simulation
     while t<tmax
   ss=ss+1;
  .............
  .................
   ................
    if ss==2000
    dlmwrite('d:\Results_for_model.txt',Tnew,'-append');
    ss=0;
    end


    end   % end of while loop

你有什么建议吗？因为在文本文件中，对于 Tnew(i,j) 值，每 10 行后，模型计算下一个 dt 值。因此，结果数据看起来像一团糟，每 10 行之后，它就会给出下一次值结果。

有没有办法根据特定的行和列来协调编写这个结果（因为否则需要一次又一次地组织大量的数据）？

我想在 3d 图中绘制这个温度分布图，在这种情况下是半球，我有 Tnew(r,theta,t)，但我对如何表示这个温度分布图以在半球图中显示感到困惑.我想听听您对此的建议。提前致谢！！

【问题讨论】：

您的循环从 9 运行到 2，因此不运行？
i = (nodes-1):2 = (J+1-1):2 = 9:2 这是一个立即终止的循环
语法是for counter = start:optional step:end
除上述之外，在代码的初始化部分中似乎没有定义节点，for i=1:nodes等
哈哈哈刚刚打败我，不错的赶上

标签： matlab transfer heat

【解决方案1】：

确切地说，MATLAB 为for loop 构造定义的正确语法是

for index = values
    program statements
    ...
end

其中values 具有以下形式之一：

initval:endval

initval:step:endval

valArray

你的代码被解析为initval:endval = 9:2，这意味着循环运行了0次，导致没有计算。

【讨论】：