[egemen]

Hi everyone,

i faced up a problem on my matlab code related to fluid flow over a horizontal plate by boundary layer theorem. I found nodal equations and wrote a script to solve it but it didnt work. I really dont understand why it did not work. I controlled it many times and i could not find any problem with code. I will be crazy. Are there someone can help me on this code?

Boundary Layer Equations for Fluid Flow

1.u(du/dx)+v(du/dy)=mu(d^2u/dy^2) (momentum eq.)
2.du/dx+dv/dy=0 (continuity eq.)

Code at below;

clear all
clc
x1=0; %initial x value
x2=0.1; %final x value (lenght of the domain)
y1=0; %initial y value
y2=0.1; %final y value
dv=1.307*10^-6; %kinematic viscosity
n=100; %number of volumes and also nodes
g=5; %initial guess
Q=g*ones(n,n); %initial Q values for each node
U=g*ones(n+1,n+1); %initial U values at surfaces
V=g*ones(n+1,n+1); %initial V values at surfaces
Uinf=4; %U infinity m/s
maxerrq=1;
maxerrv=1;
h=(x2-x1)/n; %dx=dy=h the number of nodes
D=dv/h; %kinematic viscosity over step
k=0;

while maxerrq>10^-6 && maxerrv>10^-6
maxerrq=0;
maxerrv=0;
k=k+1;
Qold=Q;

for j=1:n
U(1,j)=Uinf;
U(j,n)=Uinf;
end

for i=1:n
V(i,1)=0;
U(i,1)=0;
end

Vold=V;
Uold=U;

for i=2:n-1
for j=2:n-1
api=(U(i+1,j)+V(i,j+1)+2*D);
aei=0;
awi=U(i,j);
ani=D;
asi=D+V(i,j);
Si=0;
Q(i,j)=(aei*Q(i+1,j)+awi*Q(i-1,j)+ani*Q(i,j+1)+asi*Q(i,j-1)+Si)/api;
end
end

% corners

% lower-left corner's coefficients & nodal eq. (B1)
ap1=U(2,1)+V(1,2)+3*D;
ae1=0;
aw1=0;
an1=D;
as1=0;
S1=U(1,1)*Uinf;
Q(1,1)=(ae1*Q(2,1)+an1*Q(1,2)+S1)/ap1; %B1
% --------------------------------------------
% upper-left corner's coef. & nodal eq. (B3)
ap3=U(2,n)+3*D;
ae3=0;
aw3=0;
an3=0;
as3=V(1,n)+D;
S3=(U(1,n)-V(1,n+1)+2*D)*Uinf;
Q(1,n)=(ae3*Q(2,n)+as3*Q(1,n-1)+S3)/ap3; %B3
% --------------------------------------------
% upper-right corner's coef. & nodal eq. (B5)
ap5=U(n+1,n)+3*D;
ae5=0;
aw5=U(n,n);
an5=0;
as5=V(n,n)+D;
S5=2*D*Uinf;
Q(n,n)=(aw5*Q(n-1,n)+as5*Q(n,n-1)+S5)/ap5;
% --------------------------------------------
% lower-right corner's coef. & nodal eq. (B7)
ap7=U(n+1,1)+V(n,2)+3*D;
ae7=0;
aw7=U(n,1);
an7=D;
as7=0;
S7=0;
Q(n,1)=(aw7*Q(n-1,1)+an7*Q(n,2)+S7)/ap7;
% --------------------------------------------

% edges

% left & right edges (B2 & B6)
for j=2:n-1
Q(1,j)=(D*Q(1,j+1)+(D+V(1,j))*Q(1,j-1))/(U(2,j)+V(1,j+1)+2*D); %for B2
Q(n,j)=(U(n,j)*Q(n-1,j)+D*Q(n,j+1)+(V(i,j)+D)*Q(n,j-1))/(U(n+1,j)+V(n,j+1)+2*D); %for B6
end
% ---------------------------------------------
% lower & upper edges (B8 & B4)
for i=2:n-1
Q(i,1)=(U(i,1)*Q(i-1,1)+D*Q(i,2))/(U(i+1,1)+V(i,2)+3*D); %for B8
Q(i,n)=((V(i,n)-2*D)*Q(i,n-1)+(U(i,n)-V(i,n+1)+2*D)*Uinf)/(U(i+1,n)+3*D); %for B4
end
% ---------------------------------------------

for j=1:n
for i=2:n
U(i,j)=(Q(i-1,j)+Q(i,j))/2; %average U velocity at surfaces
end
end

for i=1:n
for j=2:n+1
V(i,j)=U(i,j-1)-U(j,i)+V(i,j-1); %continuity eq.
end
end

for i=1:n
for j=1:n
errq=abs(Qold(i,j)-Q(i,j));
errv=abs(Vold(i,j)-V(i,j));
if errq>maxerrq
maxerrq=errq;
end
if errv>maxerrv
maxerrv=errv;
end
end
end
end