[cothiele]

I want to solve Continuity and x-momentum Equations for a boundary layer over a flat plate. It is a 0 pressure gradient flow.

My equations are diverging, and I can't seem to find why they are. I am using an under-relaxation technique to obtain convergence. Could anyone help to spot my mistake?

Using a central/forward/backward difference scheme based on where I am at on my grid.

MATLAB Code is below:

%% Parameters
rho = 1.225; % kg/m^3 - Air sea level standard day
U = 10; % m/s - freestream
mu = 1.789e-5; % kg/(m*s) - Air sea level standard day
nu = mu/rho; % m^2/s
L = 0.25; % m
H = 0.05; % m

% nx and ny defines matrix grid, dx and dy for step size
nx = 51;
ny = 11;
dx = .005; % 50 grid points in x-dir - .005m spacing
dy = .005; % 10 grid points in y-dir - .005m spacing

% Re based on length
Re_L = U*L/nu;

% Filling mesh with 0 to begin for v and 10 for all u
v = zeros(ny,nx); % All initial v are 0
u = ones(ny,nx)*U; % Filling the entire mesh with all u = 10

% Boundary Conditions
% Goes (j,i) each segment, fills mesh with known values
u(:,1) = U; % Inlet
v(:,1) = 0; % Inlet
u(1, = 0; % No Slip wall
v(1, = 0; % No Slip wall
u(ny, = U; % Free stream (top of mesh)

%% Solving
% Guess values of alpha
alphau = .05; %guess, will need to figure out best value for this
alphav = .015;
% Differences
diff_u=[];
diff_v=[];
diff_u_max=[];
diff_v_max=[];
diff_both=[];

% Tolerance and convergence for main while loop
tol = .001;
z = 0; % Starting the index for iterations

while z < 100 % This is number of iterations it performs
for i = 2:1:nx-1
for j = 2:1:ny-1
if i == nx
% use backwards
v_old = v(j,i);
%v(j,i) = ((u(j,i-1) - u(j,i))/dx)*dy + (v(j-1,i));
v(j,i) = ((-u(j,i-1) - u(j,i))/dx)*dy + (v(j-1,i));
change_v = v_old - v(j,i);
v(j,i)= v_old + (alphav * change_v);
u_old = u(j,i);
u(j,i) = u(j,i-1) + (((-v(j,i) + v(j-1,i)) / dy) * dx);
change_u = u_old - u(j,i+1);
u(j,i+1) = u_old + (alphau * change_u);
diff_u = [diff_u change_u];
diff_v = [diff_v change_v];
elseif j == ny
v_old = v(j,i);
v(j,i) = v(j-1,i);
u_old = u(j,i+1);
u(j,i+1) = U;
else
% This is central differences
u_old = u(j,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);
change_u = u_old - u(j,i+1);
u(j,i+1) = u_old + (alphau * change_u);
v_old = v(j+1,i);
v(j+1,i) = v(j-1,i+1) - dy/2/dx*(u(j,i+1)-u(j,i)+u(j-1,i+1)-u(j-1,i));
change_v = v_old - v(j+1,i);
v(j+1,i)= v_old + (alphav * change_v);
diff_u = [diff_u change_u];
diff_v = [diff_v change_v];
end
end
end
diff_U = max(abs(diff_u));
diff_V = max(abs(diff_v));
diff_u_max = [diff_u_max diff_U];
diff_v_max = [diff_v_max diff_V];
diff = [diff_U diff_V];
max_diff = max(diff);

if abs(max_diff) < tol
break
end
z = z + 1;
end

% Plotting to show convergence throughout the iterations
figure
hold on
plot(diff_u_max)
plot(diff_v_max)
hold off
xlabel('Iterations');
ylabel('Chnage in u and v');
title('Convergence of u and v Throughout Iterations');
legend('u difference', 'v difference');

[FMDenaro]

The Prandtl equations are parabolic, you need to integrate along x in a time-like fashion

[cothiele]

Originally Posted by FMDenaro

The Prandtl equations are parabolic, you need to integrate along x in a time-like fashion

Isn't that what z < 100 and the z = z + 1 doing?
It's going through in a time like process where it increases 1 each time?

[FMDenaro]

Originally Posted by cothiele

Isn't that what z < 100 and the z = z + 1 doing?
It's going through in a time like process where it increases 1 each time?

if you have u(i,j) standing for u(x,y), you compute u(i+1,j) basing on the values at i. Thus, x is the time-like direction.
Be aware of the numerical stability constraint for the parabolic equation.

[cothiele]

Originally Posted by FMDenaro

if you have u(i,j) standing for u(x,y), you compute u(i+1,j) basing on the values at i. Thus, x is the time-like direction.
Be aware of the numerical stability constraint for the parabolic equation.

I'm sorry, but I am a bit confused now.

I don't normally do anything with CFD, but this is for a class project in a fluids class.

I am satisfying the stability criteria, since my step size is .005 in both dx and dy, and then my under-relaxation technique takes care of the numerical stability by
new = old + alpha*(old - calculated).

So by going u(i+1,j) I need to change my for loop from i = 2:1:nx-1 to 2:1:? if it is to be time stepped?

Thanks
[/B]

[FMDenaro]

Originally Posted by cothiele

I'm sorry, but I am a bit confused now.

I don't normally do anything with CFD, but this is for a class project in a fluids class.

I am satisfying the stability criteria, since my step size is .005 in both dx and dy, and then my under-relaxation technique takes care of the numerical stability by
new = old + alpha*(old - calculated).

So by going u(i+1,j) I need to change my for loop from i = 2:1:nx-1 to 2:1:? if it is to be time stepped?

Thanks
[/B]

The problem is solved in an open domain, you fix the initial condition at x=0 and the BCs at y=0 and y->+Inf. You have to proceed from a position x to a position x+dx by integrating along the normal to wall direction y.

There is no theoretical ending in this process, you integrate from x=0 to the position x required by your problem where you get the Re_x value.

[cothiele]

Originally Posted by FMDenaro

The problem is solved in an open domain, you fix the initial condition at x=0 and the BCs at y=0 and y->+Inf. You have to proceed from a position x to a position x+dx by integrating along the normal to wall direction y.

There is no theoretical ending in this process, you integrate from x=0 to the position x required by your problem where you get the Re_x value.

Okay. I set my boundary conditions before I started my loop, and then I have i being the x+dx. since each step is .005 it adds .005 from the previous position of i.

I am going the entire length of the plate to look at the growth of the boundary layer. from x = 0 to x = .25m, where i is each grid point I fill with a u and v velocity value

I wish I could put a photo in of my grid to show