구면 좌표계 경계 조건 구현을위한 열전달

Question

반구에 대해 열 전달 (열전도 및 대류)을 적용하고 싶습니다. 이것은 구 좌표계에서 일시적으로 균질 한 열전달입니다. 발열은 없습니다. 반구의 경계 조건은 Tinitial = 20도 실내 온도에서 시작됩니다. 외부 환경 온도는 -22도입니다. 반구가 견고한 소재라고 상상할 수 있습니다. 또한 소재가 동결 된 후 열전도율이 변하기 때문에 비선형 모델이며 온도 프로파일이 변경됩니다.

중심 온도가 -22도에 도달 할 때까지 일정한 시간 동안이 고체의 온도 프로필을 찾고 싶습니다.

이 경우 온도는 3 개의 매개 변수 T (r, theta, t)에 따라 달라집니다. 반경, 각도 및 시간. (∂T (r, θ, t))/∂t = 1/r^2 * ∂/∂r (∂T (r, θ, t))/∂r

내가 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이되고있다. 테일 온도에 도달 할 때까지 고체의 온도 프로파일을 냉각시키고 얼려 보여야합니다. 나는 그것이 왜 그렇게 변화하는지 이유를 알 수 없었다.

이 프로그램의 BC 구현에 대한 의견을 듣고 싶습니다. 또는이 조건에 따라 어떻게 변경해야합니까? 미리 감사드립니다 !!

Answer 1

코드가 너무 길어서 완전히 이해할 수 없지만 간단한 forward Euler 구성표를 사용하고있는 것처럼 보입니다. 맞습니까? 이 경우 dt이 너무 클 경우이 방법은 숫자가 unstable이 될 수 있으므로 시간 단계를 줄이려고하십시오 이렇게하면 계산 속도가 느려질 수 있습니다 (다시 말하면). 그러나 이것은 간단한 알고리즘에 대해 지불하는 가격입니다. 불안정하지 않은 alternatives methods이 있지만 방정식 시스템을 해결해야하므로 구현하기가 훨씬 어렵습니다.

오래 전이 간단한 구성표를 사용하여 열 시뮬레이션을 수행했습니다. 안정 범위 기준이 dt < (dx)^2 * c_p * rho/(6 * k) 인 것으로 나타났습니다. 여기서는 dx이 공간 단계이고, c_p은 비열이고, rho 밀도이고 k 열전도도입니다. 구형 좌표로이 값을 변환하는 방법을 모르겠습니다. 그때 내가 배운 것은 작은 시간 단계를 선택하는 것이 었습니다. 더 큰 것은 dx으로 가능한 한 큰 값으로 나타 냈습니다. dx을 2로 줄이면 상황을 안정하게 유지하기 위해 dt을 4로 줄여야합니다. 동시에 3D 문제의 경우 요소 수는 8 팩터만큼 증가합니다. 따라서 총 시뮬레이션 시간은 1/(dx)^5으로 조정됩니다 !!!