|
So I am trying to solve the heat transfer through a pipe with ODE45. I have used ODE45 to solve for 1 unknown say x before, but never for x(i) where I is a number of nodes.
My code is below. The first one sets variables and calls the RHS of the equation.
The second one is the RHS. The for loops are trying to move to the next node to be calculated, the if commands are checking boundary conditions. I think I am on the right lines, but I'm not quite there.
Script calling ODE45:
[CODE]
global AMB_COLD AMB_HOT CANDLE_T j ang slice x nq nl nr temp T length Ri Ro n_nodes Tambient T_Kill transfer Thot i a z H R K DZ DR DA T1 T2 T3 T4 T5 T6 ro volume node_mass Cp convective conductive
x=0;
%-----VARIABLES----- All the necessary variables
length = 0.1;
Ri = 0.01;
Ro = 0.015;
nr = 3;%number of radius slices
nl = 2;%numer slices along the length of the pipe
nq = 3;%angle slices
n_nodes = nr*nl*nq;
j = 0;
AMB_COLD = 298;
AMB_HOT = 398;
CANDLE_T = 698;
H = 30;
K = 100;
DZ = length/nl;
DR = (Ro-Ri)/(nr);
DA = pi()/nq;
ang = 0;
slice = 0;
ro = 8400; %density
volume = ((pi*Ro^2)-(pi*Ri^2))*length;
node_mass = ro*volume/nq;
Cp = 3770;
convective = H*DZ*DA*Ri;
conductive = (K*DZ*DA*Ri/DR);
%-----END-----
T = ones(nr,nq,nl); %Temperature array of nodes
T(:,:,:) = 298; %filling it with initial temperatures
Tx = T(:);%column vector
Q0 = Tx'; %row vector
Q = ones(nr*nl*nq,1);
[t,Q] = ode45('steady_state_transfer',[0 5],Q0);
%plot (t, Q(:,1)) %plot graph
[/CODE]
The next one is the RHS
[CODE]
function Qdot = steady_state_transfer(t, Q)
global AMB_COLD AMB_HOT CANDLE_T j ang slice x nq nl temp nr T length Ri Ro n_nodes transfer T_Kill Tambient Thot i a z H R K DZ DR DA T1 T2 T3 T4 T5 T6 ro volume node_mass Cp convective conductive
Qdot = Q*0;
for z = 1:nl
slice = slice + z*DZ;
ang = 0;
for a = 1:nq
ang = ang + a*DA;
for i = 1:nr
R = Ri + i*DR;
j = j + 1
i
z
a
if i == 1 && z~= nl && z ~= 1
%T1 = K*(DZ*DR)*(Q(j-nr)-Q(j))/(R*DA);
T2 = 0;
%T3 = K*(DZ*DR)*(Q(j+nr)-Q(j))/R*DA;
T4 = K*(DA*DZ)*(R+DR/2)*(Q(j+1)-Q(j))/DR;
T5 = K*(R*DA*DR)*(Q(j-nr*nr)-Q(j))/DZ;
T6 = K*(R*DA*DR)*(Q(j+nr*nr)-Q(j))/DZ;
elseif z ==1 && i == 1
%T1 = K*(DZ*DR)*(Q(j+nr*(nq-1))-Q(j))/(R*DA);
T2 = 0;
%T3 = K*(DZ*DR)*(Q(j+nr)-Q(j))/R*DA;
T4 = K*(DA*DZ)*(R+DR/2)*(Q(j+1)-Q(j))/DR;
T5 = H*R*DA*DR*(CANDLE_T-Q(j));
T6 = K*(R*DA*DR)*(Q(j+nr*nq+1)-Q(j))/DZ;
elseif i == 1 && z==nl
%T1 = K*(DZ*DR)*(Q(j+(nq-1)*nr)-Q(j))/(R*DA);
T2 = 0;
%T3 = K*(DZ*DR)*(Q(j+nr)-Q(j))/R*DA;
T4 = K*(DA*DZ)*(R+DR/2)*(Q(j+1)-Q(j))/DR;
T5 = K*(R*DA*DR)*(Q(j-nr*nr)-Q(j))/DZ;
T6 = H*(R*DA*DR)*(AMB_COLD-Q(j));%K*(R*DA*DR)*(Q(j+nr*nr)-Q(j))/DZ;
elseif i == nr && z~=1 && z~=nl
%T1 = K*(DZ*DR)*(Q(j-nr)-Q(j))/(R*DA);
T2 = K*DA*DZ*(R-DR/2)*(Q(j+1)-Q(j))/DR;
%T3 = K*(DZ*DR)*(Q(j+nr)-Q(j))/R*DA;
T4 = H*DA*DZ*(R+DR/2)*(AMB_COLD-Q(j));
T5 = K*(R*DA*DR)*(Q(j-nr*nr)-Q(j))/DZ;
T6 = K*(R*DA*DR)*(Q(j+nr*nr)-Q(j))/DZ;
elseif i ==nr && z ==1
%T1 = K*(DZ*DR)*(Q(j-nr)-Q(j))/(R*DA);
T2 = K*DA*DZ*(R-DR/2)*(Q(j+1)-Q(j))/DR;
%T3 = K*(DZ*DR)*(Q(j+nr)-Q(j))/R*DA;
T4 = H*DA*DZ*(R+DR/2)*(AMB_HOT-Q(j));
T5 = H*R*DA*DR*(CANDLE_T-Q(j));
T6 = K*(R*DA*DR)*(Q(j+nr*nr)-Q(j))/DZ;
elseif z == nl && i == nr
%T1 = K*(DZ*DR)*(Q(j-nr)-Q(j))/(R*DA);
T2 = K*(DA*DZ)*(R+DR/2)*(Q(j-1)-Q(j))/DR;
%T3 = K*(DZ*DR)*(Q(j+nr)-Q(j))/R*DA;
T4 = H*DA*DZ*(R+DR/2)*(AMB_COLD-Q(j));
T5 = K*(R*DA*DR)*(Q(j-nr*nr)-Q(j))/DZ;
T6 = H*(R*DA*DR)*(AMB_COLD-Q(j));
elseif z == 1 && i ~=1 && i ~= nr
%T1 = K*(DZ*DR)*(Q(j-nr)-Q(j))/(R*DA);
T2 = K*(DA*DZ)*(R+DR/2)*(Q(j-1)-Q(j))/DR;
%T3 = K*(DZ*DR)*(Q(j+nr)-Q(j))/R*DA;
T4 = K*DA*DZ*(R-DR/2)*(Q(j+1)-Q(j))/DR;
T5 = H*R*DA*DR*(CANDLE_T-Q(j));
T6 = K*(R*DA*DR)*(Q(j+nr*nr)-Q(j))/DZ;
elseif z == nl && i ~= nr && i ~= 1
%T1 = K*(DZ*DR)*(Q(j-nr)-Q(j))/(R*DA);
T2 = H*DA*DZ*(R+DR/2)*(AMB_COLD-Q(j));
%T3 = K*(DZ*DR)*(Q(j+nr)-Q(j))/R*DA;
T4 = K*DA*DZ*(R-DR/2)*(Q(j+1)-Q(j))/DR;
T5 = K*(R*DA*DR)*(Q(j-nr*nr)-Q(j))/DZ;
T6 = H*R*DA*DR*(AMB_COLD-Q(j));
else
%T1 = K*(DZ*DR)*(Q(j-nr)-Q(j))/(R*DA);
T2 = K*(DA*DZ)*(R+DR/2)*(Q(j-1)-Q(j))/DR;
%T3 = K*(DZ*DR)*(Q(j+nr)-Q(j))/R*DA;
T4 = K*DA*DZ*(R-DR/2)*(Q(j+1)-Q(j))/DR;
T5 = K*(R*DA*DR)*(Q(j-nr*nr)-Q(j))/DZ;
T6 = K*(R*DA*DR)*(Q(j+nr*nr)-Q(j))/DZ;
end
if a == 1
T1 = K*(DZ*DR)*(Q(j+nr*(nq-1))-Q(j))/(R*DA);
elseif a == nq
T3 = K*(DZ*DR)*(Q(j-nr*(nq-1))-Q(j))/R*DA;
else
T1 = K*(DZ*DR)*(Q(j-nr)-Q(j))/(R*DA);
T3 = K*(DZ*DR)*(Q(j+nr)-Q(j))/R*DA;
end
Qdot(j)= (T1 + T2 + T3 + T4 + T5 + T6)/(node_mass*Cp)
end
end
end
[/CODE]
|
Feed for this Thread