Thursday, December 14, 2006

一維環形斷面溫度分佈問題

檢視環形元素

環形底粗,剖面為梯形









梯形的各面熱傳情形

穩態情形下,熱累積率零 δT/δt=0
傳進熱流量等於傳出熱流量








環形由內而外視作等溫面擴散
而環形前後側有熱流


穩態時


熱累積率零 δT/δt=0


上述單位元素內,其熱對流熱傳導的總和即為熱累積量=0
所以由此可以建立下列平衡式:k{[t(i-1)+t(i)]/2 }{0.5*[R(i-1)+R(i)]}*(2π)}[T(i-1)–T(i)] / Δr+k{[t(I+1)+t(i)]/2 }{0.5*[R(i+1)+R(i)}*(2π)}[T(i+1)–T(i)] / Δr+2h(i)*2πr(i)Δr[Ta-T(i)] = 0


簡化


Ka=k*{[t(i-1)+t(i)]/2 };Kb=k*{[t(i+1)+t(i)]/2 }
原式可寫成
Ka[R(i-1)+R(i)]/Δr*[T(i-1)–T(i)]+Kb[R(i+1)+R(i)]/Δr*[T(i+1)–T(i)] +4h(i)r(i)Δr[Ta-T(i)]=0
再設 p(i)=Ka*[R(i-1)+R(i)]/Δr;
q(i)=Kb*[R(i+1)+R(i)]/Δr; s(i)=4*h(i)*R(i)*Δr
vpT(i-1)-(p+q+s)T(i)+qT(i+1)=-sTa i=2,3,…M-1



邊界條件


最內圈i=1,T(1)=To
最外圈i=M,多了最外環的熱流面


k{[t(M-1)+t(M)]/2 }{[R(M-1)+R(M)}/2*(2π)}[T(M-1)–T(M)] /Δr+
h(M)(2πR(M))*t(M)*[Ta-T(M)]+
2h(M)*{[R(M-1)+R(M)]/2+R(M)}/2*(2πΔr/2)*[Ta-T(M)] = 0

整理後


ss=p(M)+2h(M)R(M)t(M)+2h(M)[0.25R(M-1)+0.75R(M)]ΔrpT(M-1)-ssT(M) = -[ss-p(M)]*Ta


熱流量Qfin
R=Ri(最內圈)時,Qfin= - k(1) t(1)(2πR(1)[dT/dr]
因此從第一點流入第二點的熱流量為
Q=k(1) t(1)(2πR(1)[T1-T2]/Δr
熱對流時,與外界接觸的點溫度
Tm=[T1+Tr]/2=[T1+(T1+T2)/2]/2= (3/4)T1+(1/4)T2
Qfin=k(1) t(1)(2πR(1)[T(1)-T(2)]/Δr+ 2 (h(1)Δr/2){[(3/4)T(1) + (1/4)T(2)]-Ta}= k*t(1)*[T(1)-T(2)]/Δr + h(1)*[0.75T(1)+0.25T(2)-Ta]Δr

前置處理


n=hh(1);cc=hh(2);
mR=R/100;% length in cm
R1=mR(1);R2=mR(2);
thick=thick/1000;%
tt=linspace(thick(1),thick(2),M);%thickness at nodes
delR=(R2-R1)/(M-1);%length of each node
R=linspace(R1,R2,M);


A係數矩陣的建立
p=zeros(1,M);
q=p;s=p;T=p;T(1:M)=deal(To);
h(1:M)=deal(cc*To^0.25);(or n) h = 1.9074 (T- Ta)0.25
for i=2:M,p(i)=k*(tt(i-1)+tt(i))*(R(i-1)+R(i))/delR;endfor i=2:M-1,q(i)=k*(tt(i+1)+tt(i))*(R(i+1)+R(i))/delR;endT=T‘;h=h’; 行轉置成列

while 1h=cc*(T-Ta).^n;Tm=T(M);s=4*h.*R‘*delR;ss=p(M)+2*h(M)*R(M)*tt(M)+2*h(M)*(0.25*R(M-1)+0.75*R(M))*delR;C=[To -s(1:M-2)’*Ta -(ss-p(M))*Ta]‘;B=[0 0, q(2:M-1);1 -p(2:M-1)-q(2:M-1)-s(2:M-1)’ -ss;...p(2:M) 0]‘;A=spdiags(B,[1 0 -1],M,M);
帶狀對角矩陣以B為對角元素,[1 0 -1]為設定所在對角線位置
T=A\C;if abs(T(M)-Tm)<0.01,>使用FEMLAB模擬



執行結果與FEMLAB比較


Thursday, November 30, 2006

熱傳所有穩態邊界的結果

溫度分佈矩陣

MESH圖

SURF圖

編號

熱傳模組的建立(續) 穩態邊界所有情形


case 2
disp('Now we will divide the 2-D plane for FDE')
%too choose the divide element by setting the side
r=input('Please input m row of the matrix you want ==> ');
c=input('Please input n column of the matrix you want ==> ');
zeros(r,c);
for t=1:r*c
v(t)=t;
end
V=reshape(v,r,c);
%set up a matrix to number every element of the matrix
N=r*c;
A=zeros(N);
%after observing the FDE matrix, derive down rule code
for i=1:N
A(i,i)=-4;
if i+r0
A(i,i-r)=1;
end
if i+1> ')
H=input('The convection coeffcient is>> ')
K=input('The conduction coeffcient is>> ')
L=input('The length of element>> ')
%First of all,we deal the boundary exclude the corner
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Top boundary
switch Bu
case 1
Tu=input('Temperature of top boundary>> ')
for j=2:c-1 %it may be a bug here
B((j-1)*r+1,1)=-Tu;
end
case 2
Tu=input('The surrounding temperature of the top>> ')

for j=2:c-1
B((j-1)*r+1,1)=-(2*H*Tu*L/K+G*L^2/K);
A((j-1)*r+1,(j-1)*r+1)=-(4+2*H*L/K);% retune the matrix
A((j-1)*r+1,(j-1)*r+2)=2;
end
case 3
Qu=input('-dT/dy (toward the center is positive) =')
for j=2:c-1
B((j-1)*r+1,1)=-(2*Qu*L/K+G*L^2/K);
A((j-1)*r+1,(j-1)*r+1)=-4 ;% retune the matrix
A((j-1)*r+1,(j-1)*r+2)=2;
end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%left boundary
switch Bl
case 1
Tl=input('Temperature of left boundary>> ')
for k=2:r-1
B(k,1)=-Tl;
end
case 2
Tl=input('The surrounding temperature of the left>> ')
for k=2:r-1
B(k,1)=-(2*H*Tl*L/K+G*L^2/K);
A(k,k)=-(4+2*H*L/K);
A(k,k+r)=2;
end

case 3
Ql=input('-dT/dy (toward the center is positive) =')
for k=2:r-1
B(k,1)=-(2*Ql*L+G*L^2);
A(k,k)=-4;
A(k,k+r)=2;
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%down boundary

switch Bd
case 1
Td=input('Temperature of down boundary>> ')
for h=2:c-1
B(r*h,1)=-Td;
end
case 2
Td=input('The surrounding temperature of the bottom>> ')
for h=2:c-1
B(r*h,1)=-(2*H*Td*L/K+G*L^2/K);% retune the matrix
A(r*h,r*h)=-(4+2*H*L/K);
A(r*h,r*h-1)=2;
end
case 3
Qd=input('-dT/dy (toward the center is positive) =')
for h=2:c-1
B(r*h,1)=-(2*Qd*L/K+G*L^2/K);% retune the matrix
A(r*h,r*h)=-4;
A(r*h,r*h-1)=2;
end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%right boundary
switch Br
case 1
Tr=input('Temperature of right boundary>> ')
for p=2:r-1
B(r*(c-1)+p,1)=-Tr;
end
case 2
Tr=input('The surrounding temperature of the right>> ')
for p=2:r-1
B(r*(c-1)+p,1)=-(2*H*Tr*L/K+G*L^2/K);
A(r*(c-1)+p,r*(c-1)+p)=-(4+2*H*L/K);
A(r*(c-1)+p,r*(c-2)+p)=2;
end

case 3
Qr=input('-dT/dy (toward the center is positive) =')
for p=2:r-1
B(r*(c-1)+p,1)=-(2*H*Qr*L/K+G*L^2/K);
A(r*(c-1)+p,r*(c-1)+p)=-4;
A(r*(c-1)+p,r*(c-2)+p)=2;
end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%four corner's temp need another expression (6 situation)

% the first angle B(1,1)
if Bu==1&Bl==1 %both constant Temp
B(1,1)=-(G*L^2/(2*K)+(Tu+Tl));
A(1,1)=-4;
elseif Bu==2&Bl==2 %both convection
B(1,1)=-(G*L^2/(2*K)+H*L*(Tu+Tl)/K);
A(1,1)=-(2+2*H*L/K);
elseif Bu==3&Bl==3 %both heat flux or insulated
B(1,1)=-(G*L^2/(2*K)+2*(Qu+Ql)*L/K);
A(1,1)=-2;
elseif Bu==1&Bl==2
B(1,1)=-(G*L^2/(2*K)+H*L*Tl/K+Tu);
A(1,1)=-(3+H*L);
elseif Bu==2&Bl==1
B(1,1)=-(G*L^2/(2*K)+H*L*Tu/K+Tl);
A(1,1)=-(3+H*L);
elseif Bu==3&Bl==1
B(1,1)=-(G*L^2/(2*K)+Qu*L/K+Tl);
A(1,1)=-3;
elseif Bu==1&Bl==3
B(1,1)=-(G*L^2/(2*K)+Ql*L/K+Tu);
A(1,1)=-3;
elseif Bu==2&Bl==3
B(1,1)=-(G*L^2/(2*K)+Ql*L/K+H*L*Tu/K);
A(1,1)=-(2+H*L/K);
elseif Bu==3&Bl==2
B(1,1)=-(G*L^2/(2*K)+Qu*L/K+H*L*Tl/K);
A(1,1)=-(2+H*L/K);
end

% the second angle B(r,1)

if Bd==1&Bl==1 %both constant Temp
B(r,1)=-(G*L^2/(2*K)+(Td+Tl));
A(r,r)=-4;
elseif Bd==2&Bl==2 %both convection
B(r,1)=-(G*L^2/(2*K)+H*L*(Td+Tl)/K);
A(r,r)=-(2+2*H*L/K);
elseif Bd==3&Bl==3 %both heat flux or insulated
B(r,1)=-(G*L^2/(2*K)+2*(Qd+Ql)*L/K);
A(r,r)=-2;
elseif Bd==1&Bl==2
B(r,1)=-(G*L^2/(2*K)+H*L*Tl/K+Td);
A(r,r)=-(3+H*L);
elseif Bd==2&Bl==1
B(r,1)=-(G*L^2/(2*K)+H*L*Td/K+Tl);
A(r,r)=-(3+H*L);
elseif Bd==3&Bl==1
B(r,1)=-(G*L^2/(2*K)+Qd*L/K+Tl);
A(r,r)=-3;
elseif Bd==1&Bl==3
B(r,1)=-(G*L^2/(2*K)+Ql*L/K+Td);
A(r,r)=-3;
elseif Bd==2&Bl==3
B(r,1)=-(G*L^2/(2*K)+Ql*L/K+H*L*Td/K);
A(r,r)=-(2+H*L/K);
elseif Bd==3&Bl==2
B(r,1)=-(G*L^2/(2*K)+Qd*L/K+H*L*Tl/K);
A(r,r)=-(2+H*L/K);
end

% the third angle B(r*(c-1)+1,1)
if Bu==1&Br==1 %both constant Temp
B(r*(c-1)+1,1)=-(G*L^2/(2*K)+(Tu+Tr));
A(r*(c-1)+1,r*(c-1)+1)=-4;
elseif Bu==2&Br==2 %both convection
B(r*(c-1)+1,1)=-(G*L^2/(2*K)+H*L*(Tu+Tr)/K);
A(r*(c-1)+1,r*(c-1)+1)=-(2+2*H*L/K);
elseif Bu==3&Br==3 %both heat flux or insulated
B(r*(c-1)+1,1)=-(G*L^2/(2*K)+2*(Qu+Qr)*L/K);
A(r*(c-1)+1,r*(c-1)+1)=-2;
elseif Bu==1&Br==2
B(r*(c-1)+1,1)=-(G*L^2/(2*K)+H*L*Tr/K+Tu);
A(r*(c-1)+1,r*(c-1)+1)=-(3+H*L);
elseif Bu==2&Br==1
B(r*(c-1)+1,1)=-(G*L^2/(2*K)+H*L*Tu/K+Tr);
A(r*(c-1)+1,r*(c-1)+1)=-(3+H*L);
elseif Bu==3&Br==1
B(r*(c-1)+1,1)=-(G*L^2/(2*K)+Qu*L/K+Tr);
A(r*(c-1)+1,r*(c-1)+1)=-3;
elseif Bu==1&Br==3
B(r*(c-1)+1,1)=-(G*L^2/(2*K)+Qr*L/K+Tu);
A(r*(c-1)+1,r*(c-1)+1)=-3;
elseif Bu==2&Br==3
B(r*(c-1)+1,1)=-(G*L^2/(2*K)+Qr*L/K+H*L*Tu/K);
A(r*(c-1)+1,r*(c-1)+1)=-(2+H*L/K);
elseif Bu==3&Br==2
B(r*(c-1)+1,1)=-(G*L^2/(2*K)+Qu*L/K+H*L*Tr/K);
A(r*(c-1)+1,r*(c-1)+1)=-(2+H*L/K);
end

% the fourth angle B(r*c,1)
if Bd==1&Br==1 %both constant Temp
B(r*c,1)=-(G*L^2/(2*K)+(Td+Tr));
A(r*c,r*c)=-4;
elseif Bd==2&Br==2 %both convection
B(r*c,1)=-(G*L^2/(2*K)+H*L*(Td+Tr)/K);
A(r*c,r*c)=-(2+2*H*L/K);
elseif Bd==3&Br==3 %both heat flux or insulated
B(r*c,1)=-(G*L^2/(2*K)+2*(Qd+Qr)*L/K);
A(r*c,r*c)=-2;
elseif Bd==1&Br==2
B(r*c,1)=-(G*L^2/(2*K)+H*L*Tr/K+Td);
A(r*c,r*c)=-(3+H*L);
elseif Bd==2&Br==1
B(r*c,1)=-(G*L^2/(2*K)+H*L*Td/K+Tr);
A(r*c,r*c)=-(3+H*L);
elseif Bd==3&Br==1
B(r*c,1)=-(G*L^2/(2*K)+Qd*L/K+Tr);
A(r*c,r*c)=-3;
elseif Bd==1&Br==3
B(r*c,1)=-(G*L^2/(2*K)+Qr*L/K+Td);
A(r*c,r*c)=-3;
elseif Bd==2&Br==3
B(r*c,1)=-(G*L^2/(2*K)+Qr*L/K+H*L*Td/K);
A(r*c,r*c)=-(2+H*L/K);
elseif Bd==3&Br==2
B(r*c,1)=-(G*L^2/(2*K)+Qd*L/K+H*L*Tr/K);
A(r*c,r*c)=-(2+H*L/K);
end

%calculate and plot the distribution of temp
X=inv(A)*B;
Q=reshape(transpose(X),r,c);
figure(1)
mesh(Q);
xlabel('Top boundary')
ylabel('Left boundary')
zlabel('Temperature(℃)')
figure(2)
surf(Q);
xlabel('Top boundary')
ylabel('Left boundary')
zlabel('Temperature(℃)')

Thursday, November 23, 2006

熱傳模組的建立

解數值熱傳,所使用的是有限差分法,其方法如下:


需注意到一件事情,解數值熱傳時矩陣的編號是

有影響次序性的,因此在code中有個V矩陣是作

為編碼用的。


程式碼如下所示:


clf

clear all

format long

M=menu('請選擇熱傳的模式種類',...

'穩態無熱生成,無微分邊界',...

'穩態有熱生成,且有微分邊界',...

'非穩態且適應各種條件');

switch M

case 1

disp('Now we will divide the 2-D plane for FDE')

%too choose the divide element by setting the side

r=input('Please input m row of the matrix you want ==> ');

c=input('Please input n column of the matrix you want ==> ');

zeros(r,c);

for t=1:r*c

v(t)=t;

end

V=reshape(v,r,c);

%set up a matrix to number every element of the matrix

N=r*c;

A=zeros(N);

%after observing the FDE matrix, derive down rule code

for i=1:N

A(i,i)=-4;

if i+r0

A(i,i-r)=1;

end

if i+1> ');

for j=1:c

B((j-1)*r+1,1)=-Tu;

end



%left boundary

Tl=input('Temperature of left boundary>> ');

for k=1:r

B(k,1)=-Tl;

end



%down boundary

Td=input('Temperature of down boundary>> ');

for h=1:c

B(r*h,1)=-Td;

end



%right boundary

Tr=input('Temperature of right boundary>> ');

for p=1:r

B(r*(c-1)+p,1)=-Tr;

end



%four corner's temp need another expression

B(1,1)=-(Tu+Tl);

B(r,1)=-(Td+Tl);

B(r*(c-1)+1,1)=-(Tu+Tr);

B(r*c,1)=-(Td+Tr);



%calculate and plot the distribution of temp

X=inv(A)*B;

Q=reshape(transpose(X),r,c);

mesh(Q);

xlabel('Top boundary')

ylabel('Left boundary')

zlabel('Temperature(℃)')

case 2

disp('施工中')

otherwise

disp('施工中')

end







結果


嘗試著使用有限差分去解數值熱傳,目前來說可作穩態熱傳限於無微分邊界以及熱生成











目前只做到穩態熱傳限於無微分邊界以及熱生成

下圖為作8*8的大小時,給四邊任意溫度的熱傳分佈圖


Tuesday, June 06, 2006

小流星部落格

小流星部落格
唉 被老師抓到在後面調皮...orz 阿怎都沒人來XD

推薦一下韓劇 春天華爾滋 音樂風景都很美 劇中也有人用德文說話 應該八@@?

在奧地利取景 好美 老師 帶張風景明信片來嚐嚐如何><