|
[Sponsors] |
May 5, 2020, 11:29 |
Flat Plate Boundary Layer
|
#1 |
New Member
Colin
Join Date: Apr 2020
Posts: 7
Rep Power: 6 |
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, 12:08 |
|
#2 |
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,893
Rep Power: 73 |
The Prandtl equations are parabolic, you need to integrate along x in a time-like fashion
|
|
May 5, 2020, 12:21 |
|
#3 |
New Member
Colin
Join Date: Apr 2020
Posts: 7
Rep Power: 6 |
||
May 5, 2020, 12:28 |
|
#4 | |
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,893
Rep Power: 73 |
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 time-like direction. Be aware of the numerical stability constraint for the parabolic equation. |
||
May 5, 2020, 12:44 |
|
#5 | |
New Member
Colin
Join Date: Apr 2020
Posts: 7
Rep Power: 6 |
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 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, 13:01 |
|
#6 | |
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,893
Rep Power: 73 |
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. |
||
May 5, 2020, 13:10 |
|
#7 | |
New Member
Colin
Join Date: Apr 2020
Posts: 7
Rep Power: 6 |
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 |
||
Tags |
cfd, convergence, flat plate |
Thread Tools | Search this Thread |
Display Modes | |
|
|
Similar Threads | ||||
Thread | Thread Starter | Forum | Replies | Last Post |
LES outflow condition for a subsonic flat plate boundary layer | chris x | Main CFD Forum | 3 | August 25, 2018 06:09 |
Compressible Flat Plate with coupled hydrodynamic/thermal boundary layer | Obad | OpenFOAM Pre-Processing | 0 | September 14, 2016 15:45 |
Wrong flow in ratating domain problem | Sanyo | CFX | 17 | August 15, 2015 07:20 |
Bump in flat plate boundary layer | DCK | STAR-CCM+ | 2 | July 23, 2015 02:55 |
Flat plate boundary layer problem | student | Main CFD Forum | 3 | May 21, 2007 14:10 |