Boundary layer over a flat plate using CFD
1 Attachment(s)
I want to solve Continuity and x-momentum Equations using Finite Volume method and first 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+1Uold=U;vold=v;for i=2:cv for j=2:cvu(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; endend U(1,1)=(u(1,1)*uinf+D*U(2,1))/(u(1,2)+3*D+v(2,1)); U(cv,1)=((u(1,1)+2*D)*uinf+((v(cv,1)+D)*U(cv-1,1)))/(u(cv,2)+3*D); U(cv,cv)=(u(cv,cv)*U(cv-1,cv)+D*U(cv,cv-1)+2*D*uinf)/(u(cv,cv)+3*D); U(1,cv)=((u(1,cv))*U(1,cv-1)+D*U(2,cv))/(u(1,cv+1)+v(2,cv)+3*D); 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); 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); 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); U(i,1)=((u(i,1))*U(i-1,1)+D*U(i+1,1))/((u(i-1,1))+(v(i,2))+3*D); 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); 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 endend |
Quote:
|
no :) there is not missing line. the code looks like that when it is uploaded here, that's why I also added code.txt file..
|
Quote:
|
Just noticed this is very old thread please ignore. (You are solving 1D NS equation in 2D mesh without temperature and density??. Please discretize the equations in paper and verify.)
|
Yes - solving 1D flow with a mesh in x and y. No temperature
rho = 1.225; % kg/m^3 U = 10; % m/s - freestream mu = 1.789e-5; For a central difference I have u(i+1) = u(j,i+1) = (2*dx/u(j,i)) * (((mu/rho)*((u(j-1,i) + 2*u(j,i) + u(j+1,i))/(dy^2))) - (v(j,i)*((u(j-1,i) - u(j+1,i))/(2*dy)))) + u(j,i-1); and v(j-1) = v(j-1,i) = 2*dy*((u(j,i+1) - u(j,i-1))/(2*dx)) + v(j+1,i); |
All times are GMT -4. The time now is 07:52. |