AME535 AME 535a - Final Project Solution

Consider a steady, viscous entrance flow into a two-dimensional plane channel (see the attached figure). The channel walls are separated by a distance 2h in the y direction and at the entrance x = 0 the velocity is assumed to have a uniform value u0 and is directed in the x direction. The flow is driven by a constant pressure gradient in the streamwise direction x. Between the inflow and the outflow at x = xmax the velocity changes from the uniform value to a form appropriate for a fully developed 2-D parabolic Poiseuille flow for which
Entrance Flow

1. Using a half channel width h as a length scale, u0 as a velocity scale, and h/u0 as a time scale show that the nondimensional, 2-D Navier-Stokes equations are:
,
∂v ∂v ∂v ∂P 1 ∂2v ∂2v
= −u − v − + 2 + ∂y2!,
∂t ∂x ∂y ∂y Re ∂x
∂u ∂v
+ = 0, ∂x ∂y
where u = u(x,y) and v = v(x,y) are the velocity components, P(x,y) is the pressure, and Re = u0h/ν is the Reynolds number (ν is the kinematic viscosity).
2. The pressure is decomposed into a mean pressure P¯(x) driving the flow and a perturbation p(x,y), i.e. P = P¯(x)+p(x,y). Show that at the outflow the velocity field is
,
where |C| is the constant mean pressure gradient, i.e . Note that the system of coordinates is chosen such that −1 < y < +1.
3. Use the conservation of mass to show that the constant |C| = 3/Re and that the u-momentum equation can be re-written as
.
4. Write a numerical code to determine the flow in the channel assuming no-slip boundary conditions on the channel walls: u(x,−1) = 0, u(x,+1) = 0; a uniform streamwise velocity at the inflow: u(0,y) = 1, v(0,y) = 0; a Neumann boundary condition for the streamwise velocity at the outflow = 0 and v(xmax,y) = 0.
Method of solution:
Use a pseudo-transient, fractional step method (time-splitting method) consisting of the following steps.
1. Nonlinear step.
,
,
where (Fu,Fv) are the nonlinear terms
∂u ∂u 3
Fu = −u − v + ,
∂x ∂y Re
∂v ∂v
Fv = −u − v .
∂x ∂y
Time advancement using Adams-Bashforth method, i.e.
,
.
2. Pressure step.
∂u ∂p
= − , ∂t ∂x
∂v ∂p
= − , ∂t ∂y
∂u ∂v + = 0.
∂x ∂y
Time discretization
u∗∗ − u∗
= −∇p. ∆t
Taking divergence of the above equation and assuming that u∗∗ is incompressible, i.e. = 0, we get a Poisson equation for p !
The above equation for p is first solved with the Neumann boundary conditions on all boundaries
,
and
,
and then u∗∗ is found from the previous equation
,
.
3. Two viscous steps.
The first viscous step

.
Discretization using Crank-Nicolson scheme, i.e. for u
,
and a similar equation for v∗∗∗.
The second viscous step similarly solves equations
∂u 1 ∂2u
= 2
∂t Re ∂y
∂v 1 ∂2v
= 2,
∂t Re ∂y
using the Crank-Nicolson method. For u we get
,
and a similar equation for vn+1.
AME 535a Final Project - Hints
• For spatial derivatives use the second order central difference formulas.
• Velocity equations are solved on the internal points j = 2,3,...,NX−1 and k = 2,3,...,NY − 1 and the outflow boundary j = NX.
• Maximum time step ∆t is estimated from the CFL condition using the inflow velocity.
• For pressure formulate and solve discretized equation on all mesh points j = 1,2,...,NX and k = 1,2,...,NY assuming that the r.h.s. of the equation is zero at all boundaries. To get spatial derivatives at the boundary points, when needed, fictitious points (usually denoted with a subscript ‘0’ in the class notes) outside the computational domain must be introduced.
.

• It is preferable to construct the matrix in the direct solver for the pressure equation such that the solution is obtained as a sequence of values in columns along y, rather than in rows along x, i.e.
(p1,1,p1,2,p1,3,...,pj,k,pj,k+1,pj,k+2,...,pnx,ny−1,pnx,ny)
• The viscous step solution is similar to that used in program TWDIF.
• Test each subroutine separately.
• In order to test subroutines a good approach is to consider a so-called method of manufactured solutions (artificial exact solutions). For instance, for the pressure equation the exact solution for a square domain 0 < x < 1, 0 < y < 1 can be constructed as follows. The equation is
.
and the boundary conditions are ∂x∂p = 0 and ∂y∂p = 0 at vertical and horizontal boundaries (x=0,1; y=0,1), respectively. Consider any function that would satisfy such boundary conditions, e.g.
pex(x,y) = cos(πx)cos(πy).
For this function the Laplacian can be computed

Therefore the rhs can be used as the function f(x,y) in the equation for p(x,y), the equation solved using a given numerical method, and the numerical solution compared with pex above to determine how accurate the numerics are.
• Consider again the above pressure equation for a square domain

with the boundary conditions = 0 and ∂y∂p = 0 at vertical and horizontal boundaries (x=0,1; y=0,1), respectively. If one adds any linear function of x and y to the solution p(x,y)
p0(x,y) = p(x,y) + αx + βy + γ,
the function p0(x,y) is also a solution. To enforce the boundary conditions we must have α = 0 and β = 0 but γ is an arbitrary constant. If subroutines FACT and SOLVE are used to solve the discretized pressure equation, usually the solution is dominated by a very large constant γ. However, the derivatives ∂x∂p and ∂y∂p needed to advance in time the velocity are finite and well behaved (though not very accurate).
To get the finite values of pressure the solver can be modified by forcing the pressure to a constant value at one or few selected points. For instance, for the example of the exact solution, pex(1,1) = 1. The matrix equation can be considered for all points except point (NX,NY), corresponding for a square domain to x = 1, y = 1. In all remaining equations wherever pNX,NY appears, it is replaced by its known value of 1, and moved to the r.h.s. The result is a system of (NX ×NY −1) equations for pressure values pj,k at all points except the corner point (NX,NY ). Solving this system by subroutines FACT/SOLVE produces values of pressure of the order unity. However, the pressure gradients in both methods are approximately the same, leading to comparable results for the velocities u and v.
The approach has been tested using FACT/SOLVE and Matlab direct inverse INV, producing results with the following accuracy:
for FACT()/SOLVE()
mesh: 11X11 21X21 31X31
rms : 9.4150E-03 2.3240E-03 1.02890E-03
==============================================================
for Matlab (direct solver: inv() )
mesh 11X11 21X21 31X31
rms : 9.1395E-02 4.7675E-02 3.2274E-02
Format of the Final Report
The final report must provide the following information:
• All equations for the algorithms used in the program, including discretized equations that are solved in each step of the method.
• Solutions for three different Reynolds numbers Re = 1, Re = 10, and Re = 100. Experiment with different mesh sizes and values of xmax until solution at the outflow converges to the form predicted analytically.
• For a uniform comparison between different programs set the length of the domain in the x-direction to xmax = 3, the mesh size to ∆x = ∆y = 0.1, i.e. number of mesh points is nx = 31 and ny = 21, and the time step to ∆t = 0.01. The domain size is a bit too small for the case Re = 100 so consider only two values of Re, 1 and 10 for this part of the project.
• Provide printout and plots for the following quantities:
1. vertical profiles of the velocity components u and v at j = 2,(nx+ 1)/2,nx − 1,nx, i.e., (u(2,k),k = 1,ny), etc.; plot profiles for u in one figure together with the exact solution at the outflow.
2. velocity u and v along the centerline k = (ny+1)/2, i.e, (u(j,(ny+ 1)/2),j = 1,nx), etc.; plot u as a function of x and mark in the figure the asymptotic value of 1.5 at the outflow.
3. the flow rate Q(x) in the channel at each j; plot this quantity together with a line corresponding to the constant flow rate at the inflow.
4. the rms error for the u velocity prediction at the outflow

5. the rms error for the flow rate prediction along the channel

• After completing the above tasks explore other choices of the mesh size, domain size, time step, etc. until your program, in your opinion, gives the best possible agreement with the theoretical u velocity profile at the outflow. Print and plot your best results and compare with the results obtained for previously prescribed resolution. Your best overall results for each Reynolds number should be included in the report as well.