Computational Fluid Dynamics, flat plate boundary layer

5 views (last 30 days)
[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
  4 Comments
Matt Tearle
Matt Tearle on 13 May 2011
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
Florin Neacsu on 16 May 2011
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.

Sign in to comment.

Answers (2)

Florin Neacsu
Florin Neacsu on 16 May 2011
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

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

Categories

Find more on Computational Fluid Dynamics (CFD) in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!