Computational Fluid Dynamics, flat plate boundary layer

[EDIT: Thu May 12 23:44:55 UTC 2011 Reformat - MKF]

I want to solve Continuity and x-momentum Equations using Finite Volume method and upwind scheme. I wrote a code for this problem but it doesnt work. I got NaN error. Can you help me to find my mistake? Here is my code:

clear all
clc
cv=10;
xi=0;
xf=0.1;
yi=0;
yf=0.1;
dx=(xf-xi)/cv;
dy=(yf-yi)/cv;
U=zeros(cv);
v=zeros(cv+1);
u=zeros(cv+1);
uinf=1;
maxerr=1;
k=0;
nu=15.11*10^-6;
D=nu/dy;
u(1,:)=0;

while maxerr>0.000001

    k=k+1
Uold=U;
vold=v;
for i=2:cv
    for j=2:cv
u(i,1)=uinf;
u(i,cv+1)=u(i,cv);
u(i,j)=(U(i-1,j)+U(i,j))/2;
u(cv+1,1)=uinf;
u(cv+1,i)=uinf;
u(i,cv+1)=u(i,cv);
u(1,cv+1)=u(1,cv);
u(i,cv+1)=u(i,cv);
u(cv+1,cv+1)=u(cv+1,cv);
u(1,j)=0;
    end
end
        U(1,1)=(u(1,1)*uinf+D*U(2,1))/(u(1,2)+3*D+v(2,1)); %boundary7
        U(cv,1)=((u(1,1)+2*D)*uinf+((v(cv,1)+D)*U(cv-1,1)))/(u(cv,2)+3*D); %boundary8
        U(cv,cv)=(u(cv,cv)*U(cv-1,cv)+D*U(cv,cv-1)+2*D*uinf)/(u(cv,cv)+3*D); %boundary4
        U(1,cv)=((u(1,cv))*U(1,cv-1)+D*U(2,cv))/(u(1,cv+1)+v(2,cv)+3*D);%boundary3        

    for i=2:cv-1
            for j=2:cv-1

                U(1,j)=(((v(1,j-1))+D)*U(i-1,j)+u(1,1)*uinf)/(u(i,j)+(v(i,j))+2*D);%boundary2
                U(cv,j)=((u(cv,j))*U(cv-1,j)+((v(cv,j))+D)*U(cv,j-1)+D*U(cv-1,j))/(u(cv,j)+(v(cv,j))+2*D); %boundary5
                U(i,cv)=(((u(i,cv)))*U(i-1,cv)+((v(i,cv))+D)*U(i,cv-1)+2*D*uinf)/(u(i,cv+1)+3*D); %boundary9
                U(i,1)=((u(i,1))*U(i-1,1)+D*U(i+1,1))/((u(i-1,1))+(v(i,2))+3*D); %boundary6
                U(i,j)=(u(i-1,j)*U(i-1,j)+(v(i,j-1)+D)*U(i,j-1)+D*U(i,j+1))/(u(i+1,j)+v(i,j+1)+2*D); %internalnodes
        end
    end
    for i=2:cv
        for j=2:cv-1
            v(i,j+1)=u(i+1,j)-u(i-1,j)+u(i,j-1);
        end
    end
     for i=1:cv
          for j=1:cv
              er=abs(U(i,j)-Uold(i,j));              
              if er>maxerr
                 maxerr=er;
              end
          end         
     end
end

[Florin Neacsu]

Hi

Try this :

er=abs(U(i,j)-Uold(i,j));

if isnan(er)

pause

end

if er>maxerr

maxerr=er;

end

And insert a break point on the line that says "pause". You will see that you U and Uold have value in the range of 1e160 and 1e198. So I guess your eq. is diverging. Check your code again.

And from a logical point of view : (correct me if I'm wrong)

you start with maxer=1

your big loop is maxerr>0.000001

if er>maxerr

maxerr=er

So you only change maxerr if er>maxerr. But if maxerr=1 that means you only change maxerr to something grater than 1, so your while loop never ends ?

Regards,

Florin

[Andrew Newell]

Why not make this an answer and collect some votes, Florin?

[Matt Tearle]

Also so hs can accept your answer and we can all rest easy (b/c it looks like you're all over this).

[Florin Neacsu]

Hi,

Indeed I should have posted as an answer. But, OP hasn't replied in 4 days so I figure he/she either fixed his/her code or he/she's not that interested in it.

Answer 1

[Florin Neacsu]

Hi

Try this:

er=abs(U(i,j)-Uold(i,j));
if isnan(er)
pause
end
if er>maxerr
maxerr=er;
end

And insert a break point on the line that says "pause". You will see that you U and Uold have value in the range of 1e160 and 1e198. So I guess your eq. is diverging. Check your code again.

And from a logical point of view : (correct me if I'm wrong)

you start with maxer=1

your big loop is maxerr>0.000001

if er>maxerr
maxerr=er

So you only change maxerr if er>maxerr. But if maxerr=1 that means you only change maxerr to something grater than 1, so your while loop never ends ?

Regards, Florin

Answer 2

[hs]

Wow, thank you a lot Florin. I couldn't solve my problem yet but I'll definitely try your opinion also.