1

半球に熱伝達（熱伝導と対流）を適用したい。これは、球座標における過渡的な均一熱伝達です。発熱はありません。半球の境界条件は、Tinitial=20℃の室温で始まります。外気温はマイナス30度。半球はソリッドな素材であることが想像できます。また、材料が凍結した後に熱伝導率が変化し、温度プロファイルが変化するため、これは非線形モデルです。

中心温度が-30℃に達するまでの特定の時間中のこの固体の温度プロファイルを見つけたいです。

この場合、Temperature は 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;

EDIT *おかげで、スクリプトが機能し、計算されるようになりました。そして、モデル システムの温度プロファイルを確認できます。

ただし、この半球の温度プロファイルを 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) がありますが、この温度プロファイルを半球グラフに表示する方法に混乱しています。それについてあなたの提案を聞きたいです。前もって感謝します ！！

4

1 に答える 1

1

正確に言うと、MATLAB で定義されている構造体の適切な構文は次のとおりfor loopです。

for index = values
    program statements
    ...
end

wherevaluesは、次のいずれかの形式です。

initval:endval

initval:step:endval

valArray

コードはinitval:endval=9:2に解析されます。これは、ループが 0 回実行され、計算が行われないことを意味します。

于 2013-04-09T15:35:40.257 に答える