Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
ODE45 with arrays

Subject: ODE45 with arrays

From: Alex

Date: 5 Apr, 2010 17:55:27

Message: 1 of 1

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]

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us