Flat Plate Boundary Layer
I want to solve Continuity and xmomentum 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 underrelaxation 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.789e5; % 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 xdir  .005m spacing dy = .005; % 10 grid points in ydir  .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:nx1 for j = 2:1:ny1 if i == nx % use backwards v_old = v(j,i); %v(j,i) = ((u(j,i1)  u(j,i))/dx)*dy + (v(j1,i)); v(j,i) = ((u(j,i1)  u(j,i))/dx)*dy + (v(j1,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,i1) + (((v(j,i) + v(j1,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(j1,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(j1,i) + 2*u(j,i) + u(j+1,i))/(dy^2)))  (v(j,i)*((u(j1,i)  u(j+1,i))/(2*dy)))) + u(j,i1); 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(j1,i+1)  dy/2/dx*(u(j,i+1)u(j,i)+u(j1,i+1)u(j1,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'); 
The Prandtl equations are parabolic, you need to integrate along x in a timelike fashion

Quote:
It's going through in a time like process where it increases 1 each time? 
Quote:
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 timelike direction. Be aware of the numerical stability constraint for the parabolic equation. 
Quote:
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 underrelaxation 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:nx1 to 2:1:? if it is to be time stepped? Thanks [/B] 
Quote:
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. 
Quote:
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 
All times are GMT 4. The time now is 09:17. 