CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Main CFD Forum

Need Help in Debugging Matlab Code

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Display Modes
Old   August 14, 2009, 12:05
Post Need Help in Debugging Matlab Code
  #1
New Member
 
Prahallada
Join Date: Aug 2009
Posts: 2
Rep Power: 0
Prahallada is on a distinguished road
Hi guys...Im trying to formulate a mathematical model, which was already proposed in a paper using matlab for my B.Tech final year project work.I finished most of the coding part and when i try to run the code i get an error,which i'm unaware of.

Code:
clear
%chan = ddeinit('excel', 'Sheet1');
R=.02%input('enter the value of  inside radius :');
Ro=.03;
l=0.5%input('enter the value of length of heat source:');
%input('enter the length of pipe:');
l1=0.1%input('enter lengh between two heat source:');
Te=20%input('enter the value of free stream temperature:');
q=1000%input('enter the value of heat flux:');
U=5%input('enter free stream velocty:');
k=12%input('enter the value of thermal conductivity ofi inside fluid:');
kw=12;
alpha=.01%input('enter the value of thermal difusivity of inside fluid:');
alpha1=4*alpha;
a=input('enter no of grid in radial directon must be odd no:');
b=input('enter no of grid in axial direction :');
n=input('enter no of heat sources on periphery:');
er=input('enter value of max possible percentage error:');
s=input('enter relaxatation factor around 1.05:')
Pe=U*2*R/alpha
deln1=(Ro/R)/(a-1);
L=n*(l+l1);
c=L/(R*(b-1));
K=kw/k;
A1=alpha/alpha1;
h=0;
delt(1,1)=0;
Nu=0;
g=((a-1)*(R/Ro)+1)
warning off all
  for i=1:a
             for j=1:b
                
                % T(i,j)=Te;
                 %t(i,j)=(T(i,j)-Te)/(q*R/k);
                 t(i,j)=0.1;
             end
   end
 format long g
 erfmx=10;
 erf=1;
 erf1=0;
 itn=1;
while(er<erfmx)
     erfmx=0;
    for i=1:a
     for j=1:b
          if j==1
          t(i,1)=0;
          end
           tn(i,b)=0;
        if j>1 && j<b
            if i==1
                X=Pe*c + 1;
                Y=Pe*c + 2;
                tn(i,j)=(t(i,j+1)+X*t(i,j-1))/Y;
            end
            if i>1 && i<((a-1)*(R/Ro)+1)
                  A=(Pe*(1-(((i-1)*Ro)/((a-1)*R))^2))/c + 2*((a-1)*(R/Ro))^2 + 2/c^2 + (((a-1)*(R/Ro))^2)/(i-1);
                  B=((a-1)*(R/Ro))^2  - (((a-1)*(R/Ro))^2)/(2*(i-1));
                  C=((a-1)*(R/Ro))^2 + (((a-1)*(R/Ro))^2)/(2*(i-1));
                  D=(Pe*(1-(((i-1)*Ro)/((a-1)*R))^2))/c + 1/c^2;
                  E=1/c^2;
             tn(i,j)=(B*t(i-1,j)+C*t(i+1,j)+D*t(i,j-1)+E*t(i,j+1))/A;
            end
            if i==((a-1)*(R/Ro)+1)
                A=(Pe*(1-(((i-1)*Ro)/((a-1)*R))^2))/c + 2*((a-1)*(R/Ro))^2 + 2/c^2 + (((a-1)*(R/Ro))^2)/(i-1);
                B=((a-1)*(R/Ro))^2  - (((a-1)*(R/Ro))^2)/(2*(i-1));
                C=((a-1)*(R/Ro))^2 + (((a-1)*(R/Ro))^2)/(2*(i-1));
                D=(Pe*(1-(((i-1)*Ro)/((a-1)*R))^2))/c + 1/c^2;
                E=1/c^2;
                x=(L/(b-1))*(j-1);
                o=3;
                if   x>(o-1)*(l+l1) && x<=(o*l+(o-1)*l1)
                 tn(i,j)=(B*t(i-1,j)+C*(t(i-1,j)+2/(a-1))+D*t(i,j-1)+E*t(i,j+1))/A; 
                end
                 if  x<=(o-1)*(l+l1)  &&  x>(o*l+(o-1)*l1) 
                   tn(i,j)=((B+C)*t(i-1,j)+D*t(i,j-1)+E*t(i,j+1))/A;
                 end  
            end 
              if i>((a-1)*(R/Ro)+1) && i<a
                  A=(Pe*(1-(((i-1)*Ro)/((a-1)*R))^2))/c + 2*((a-1)*(R/Ro))^2 + 2/c^2 + (((a-1)*(R/Ro))^2)/(i-1);
                  B=((a-1)*(R/Ro))^2  - (((a-1)*(R/Ro))^2)/(2*(i-1));
                  C=((a-1)*(R/Ro))^2 + (((a-1)*(R/Ro))^2)/(2*(i-1));
                  D=(Pe*(1-(((i-1)*Ro)/((a-1)*R))^2))/c + 1/c^2;
                  E=1/c^2;
                  tn(i,j)=(B*t(i-1,j)+C*t(i+1,j)+E*(t(i,j-1)+t(i,j+1)))/(2*((a-1)*(R/Ro))^2 + 2/c^2);
              end
                  if i==a
                    A=(Pe*(1-(((i-1)*Ro)/((a-1)*R))^2))/c + 2*((a-1)*(R/Ro))^2 + 2/c^2 + (((a-1)*(R/Ro))^2)/(i-1);
                    B=((a-1)*(R/Ro))^2  - (((a-1)*(R/Ro))^2)/(2*(i-1));
                    C=((a-1)*(R/Ro))^2 + (((a-1)*(R/Ro))^2)/(2*(i-1));
                    D=(Pe*(1-(((i-1)*Ro)/((a-1)*R))^2))/c + 1/c^2;
                    E=1/c^2;
                    x=(L/(b-1))*(j-1);
                    o=3;
                      if   x>(o-1)*(l+l1) && x<=(o*l+(o-1)*l1)
                          tn(i,j)=(B*t(i-1,j)+C*(Ro/(R*(a-1))+t(i-1,j))+E*(t(i,j-1)+t(i,j+1)))/(2*((a-1)*(R/Ro))^2 + 2/c^2);
                      end
                      if  x<=(o-1)*(l+l1)  &&  x>(o*l+(o-1)*l1) 
                      tn(i,j)=((B+C)*t(i-1,j)+E*(t(i,j-1)+t(i,j+1)))/(2*((a-1)*(R/Ro))^2 + 2/c^2);
                      end  
                  end
        end
        if j==b
             if i==1
             tn(i,j)=t(1,j-1);
             end
             if i>1 && i<((a-1)*(R/Ro)+1)
                  A=(Pe*(1-(((i-1)*Ro)/((a-1)*R))^2))/c + 2*((a-1)*(R/Ro))^2 + 2/c^2 + (((a-1)*(R/Ro))^2)/(i-1);
                  B=((a-1)*(R/Ro))^2  - (((a-1)*(R/Ro))^2)/(2*(i-1));
                  C=((a-1)*(R/Ro))^2 + (((a-1)*(R/Ro))^2)/(2*(i-1));
                  D=(Pe*(1-(((i-1)*Ro)/((a-1)*R))^2))/c + 1/c^2;
                  E=1/c^2;
                  tn(i,j)=(B*t(i-1,j)+(D+E)*t(i,j-1)+ C*t(i+1,j))/A;
             end
             if i==((a-1)*(R/Ro)+1)
                  A=(Pe*(1-(((i-1)*Ro)/((a-1)*R))^2))/c + 2*((a-1)*(R/Ro))^2 + 2/c^2 + (((a-1)*(R/Ro))^2)/(i-1);
                  B=((a-1)*(R/Ro))^2  - (((a-1)*(R/Ro))^2)/(2*(i-1));
                  C=((a-1)*(R/Ro))^2 + (((a-1)*(R/Ro))^2)/(2*(i-1));
                  D=(Pe*(1-(((i-1)*Ro)/((a-1)*R))^2))/c + 1/c^2;
                  E=1/c^2;
                  tn(i,j)=((B+C)*t(i-1,j)+(D+E)*t(i,j-1))/A;
             end
              if i>((a-1)*(R/Ro)+1) && i<a
                  A=(Pe*(1-(((i-1)*Ro)/((a-1)*R))^2))/c + 2*((a-1)*(R/Ro))^2 + 2/c^2 + (((a-1)*(R/Ro))^2)/(i-1);
                  B=((a-1)*(R/Ro))^2  - (((a-1)*(R/Ro))^2)/(2*(i-1));
                  C=((a-1)*(R/Ro))^2 + (((a-1)*(R/Ro))^2)/(2*(i-1));
                  D=(Pe*(1-(((i-1)*Ro)/((a-1)*R))^2))/c + 1/c^2;
                  E=1/c^2;
                  tn(i,j)=(B*t(i-1,j)+C*t(i+1,j)+2*E*t(i,j-1))/(2*((a-1)*(R/Ro))^2 + 2/c^2);
              end
                  if i==a
                    A=(Pe*(1-(((i-1)*Ro)/((a-1)*R))^2))/c + 2*((a-1)*(R/Ro))^2 + 2/c^2 + (((a-1)*(R/Ro))^2)/(i-1);
                    B=((a-1)*(R/Ro))^2  - (((a-1)*(R/Ro))^2)/(2*(i-1));
                    C=((a-1)*(R/Ro))^2 + (((a-1)*(R/Ro))^2)/(2*(i-1));
                    D=(Pe*(1-(((i-1)*Ro)/((a-1)*R))^2))/c + 1/c^2;
                    E=1/c^2;
                    tn(i,j)=((B+C)*t(i-1,j)+2*E*t(i,j-1))/(2*((a-1)*(R/Ro))^2 + 2/c^2);
                   
                  end
        end
        erf=abs(((tn(i,j)-t(i,j))/(tn(i,j)+t(i,j)))*200);
                t(i,j)=(1-s)*t(i,j)+s*tn(i,j);
                     if erf>erfmx
                        erfmx=erf;
                     end   
      itn=itn+1;   
     end
    end
        
        itn=itn+1;
end
itn;
t';
erfmx;
t;
T= t*(q*R/k)+Te;
 itn 
 h=(Ro/R)/(a-1);
 v(1)=0;
f(1)=0;
s=f(1);
for j=1:b
for i=2:2:g-2
    v(i)=(i-1)*(Ro/R)/(a-1);
    f(i)=(v(i)-v(i)^3)*t(i,j);
    s=s+4*f(i);
    v(i+1)=(i)*(Ro/R)/(a-1);
    f(i+1)=(v(i+1)-v(i+1)^3)*t(i+1,j);
    s=s+2*f(i+1);
end
i=g-1;
v(g-1)=(g-2)*h;
p=(g-2)*h;
s=s+4*(v(g-1)-(v(g-1)^3))*t(i,j);
s=s;
sum=1*s/(3*(g-1));
i=g;
tb(i,j)=4*sum;
delt(i,j)=t(g,j)-tb(i,j);
Nu(i,j)=0;
x=(L/(b-1))*(j-1);
             o=3;
    if x>(o-1)*(l+l1)&& x<= (o*l+(o-1)*l1)
        Q(i,j)=1;
        Nu(i,j)=2*(Ro/R)*Q(i,j)/delt(i,j);
    end
end
Nu;
tb;
delt;

fid=fopen('D:\matlabresult\hs1.txt','wt')
          fprintf(fid,'radius    nonrad   inc   axis      ndaxis      incx     temp       nondt     Nu     delt     tb\n');
for i=1:a
       n1(i)=(i-1)/(a-1);
       r(i)=n1(i)*R;
   
 for j=1:b
        x(j)= (L/(b-1))*(j-1);
          xi(j)=x(j)/R;
         rt=[r(i);n1(i);i;x(j);xi(j);j;T(i,j);t(i,j);Nu(i,j);delt(i,j);tb(i,j)];
        % csvwrite('rk.txt',rt); 
          %[c,h]=contour(x,r,T);
          %set(h,'ShowText','on','TextStep',get(h,'LevelStep')*2)
          %colormap cool
          %o=((i-1)*b)+j;
          %xlabel('x')
          %ylabel('r')
         %clabel(cs)
         
         fprintf(fid,'%1.4f   %2.4f  %3d    %1.4f     %2.4f     %3d       %6.4f    %1.4f    %1.2f    %1.2f    %1.2f\n',rt);
       
      
        
  end
 end
fclose(fid);
I'm getting an error in the line which is bolded and coloured red in the code.The error states "Subscript indices must either be real positive integers or logicals".

When i use a value of 0.02 for Ro the code is working fine without any errors.The error starts popping only for Ro=0.03.
Can some one please help me out in debugging this code.I would be thankful if some one can do so

Thanking You
Prahallada J

Last edited by Prahallada; August 14, 2009 at 12:22.
Prahallada is offline   Reply With Quote

Old   August 14, 2009, 13:27
Default
  #2
f-w
Senior Member
 
Join Date: Apr 2009
Posts: 129
Rep Power: 9
f-w is on a distinguished road
v(g-1)=(g-2)*h

In the line above, (g-1) must be a real, positive integer (it can't be 1.222 or 4.11).

Your formulations for g, R, Ro, a:
a= an integer
R=.02
Ro=.03;
g=((a-1)*(R/Ro)+1)

So, when Ro=0.02, "g" works out to be an integer, otherwise, it's not. So, make it a multiple of R.
f-w is offline   Reply With Quote

Old   August 15, 2009, 03:30
Default
  #3
New Member
 
Prahallada
Join Date: Aug 2009
Posts: 2
Rep Power: 0
Prahallada is on a distinguished road
Thanks a lot dude.....
Prahallada is offline   Reply With Quote

Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
MATLAB fractional step code Darren Main CFD Forum 6 June 16, 2013 07:21
Debugging Unsteady 2-D Panel Method Code: Wake Modeling RajeshAero Main CFD Forum 5 November 10, 2011 06:48
Matlab Code shosho Main CFD Forum 3 August 8, 2009 11:00
Running Matlab script from Scheme code dhimans Fluent UDF and Scheme Programming 0 July 28, 2009 15:13
Design Integration with CFD? John C. Chien Main CFD Forum 19 May 17, 2001 15:56


All times are GMT -4. The time now is 12:29.