Flat Plate Boundary Layer

 User Name Remember Me Password
 Register Blogs Members List Search Today's Posts Mark Forums Read

 May 5, 2020, 10:29 Flat Plate Boundary Layer #1 New Member   Colin Join Date: Apr 2020 Posts: 7 Rep Power: 5 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');

 May 5, 2020, 11:08 #2 Senior Member   Filippo Maria Denaro Join Date: Jul 2010 Posts: 6,763 Rep Power: 71 The Prandtl equations are parabolic, you need to integrate along x in a time-like fashion

May 5, 2020, 11:21
#3
New Member

Colin
Join Date: Apr 2020
Posts: 7
Rep Power: 5
Quote:
 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?

May 5, 2020, 11:28
#4
Senior Member

Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,763
Rep Power: 71
Quote:
 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.

May 5, 2020, 11:44
#5
New Member

Colin
Join Date: Apr 2020
Posts: 7
Rep Power: 5
Quote:
 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]

May 5, 2020, 12:01
#6
Senior Member

Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,763
Rep Power: 71
Quote:
 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.

May 5, 2020, 12:10
#7
New Member

Colin
Join Date: Apr 2020
Posts: 7
Rep Power: 5
Quote:
 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

 Tags cfd, convergence, flat plate

 Thread Tools Search this Thread Search this Thread: Advanced Search Display Modes Linear Mode

 Posting Rules You may not post new threads You may not post replies You may not post attachments You may not edit your posts BB code is On Smilies are On [IMG] code is On HTML code is OffTrackbacks are Off Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post chris x Main CFD Forum 3 August 25, 2018 05:09 Obad OpenFOAM Pre-Processing 0 September 14, 2016 14:45 Sanyo CFX 17 August 15, 2015 06:20 DCK STAR-CCM+ 2 July 23, 2015 01:55 student Main CFD Forum 3 May 21, 2007 13:10

All times are GMT -4. The time now is 01:36.

 Contact Us - CFD Online - Privacy Statement - Top