1

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

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

この場合、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を使って有限差分法を適用しましたが、境界条件に問題があります。半球の表面には対流があり、内側のノードでは伝導があり、半球の底は一定の温度、つまり気温 (-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 を実装するための提案、またはこの条件に従ってどのように変更する必要があるかについてお聞きしたいと思います。前もって感謝します ！！

4

1 に答える 1

1

あなたのコードは長すぎて完全には理解できませんが、単純なフォワード オイラースキームを使用しているように見えますが、正しいですか? が大きすぎるとdt、この方法は数値的に不安定になる可能性があるためです。dtこれにより、計算速度が遅くなる可能性がありますが (これも大幅に)、このような単純なアルゴリズムに対して支払う代償です。不安定性に悩まされない代替方法がありますが、連立方程式を解く必要があるため、実装がはるかに困難です。

私はずっと前に、この単純なスキームを使用していくつかの熱シミュレーションを行いました。安定性基準は であることがわかりましたdt < (dx)^2 * c_p * rho / (6 * k)。これは 3D デカルト グリッドでのシミュレーションに有効であり、dxは空間ステップであり、材料のc_p比熱、rho密度、およびk熱伝導率です。これを球座標であなたのケースに変換する方法がわかりません。そこで私が学んだことは、小さな時間ステップを選択することでしたが、より重要なdxことはできるだけ大きな時間ステップを選択することでした: 2 分の 1 に減らす場合は、物事を安定させるために 4 分の 1dxに減らす必要もあります。dt同時に、3D 問題の場合、要素の数は 8 倍になります。したがって、合計シミュレーション時間は1 / (dx)^5!!!でスケーリングされます。

于 2013-08-29T19:34:55.623 に答える