You are on page 1of 4

Final Project-Continuum Mechanics

ETH Mechanics of Continua FS 2014

November 14, 2016

(Translation Keith Patarroyo)

Exercise 10.4 Numerical solution of a stationary 2-dimensional incom-

pressible viscous flow
In this task we are looking for a steady 2-dimensional incompressible viscous flow through a straight delimited
fluid pipe, in which a rectangular object is found, see Fig. 1. Therefore, the Navier-Stokes equation to be solved
together with the continuity equation,

 + (v )v = 2 v 1 p,
 and v =0 (1)
We solve this non-trivial task hydrodynamic numerically. For this we give away the skeleton of the model
program, where only the boundary conditions and the discretized equations must be supplemented.

Figure 1: Geometry of the task, note the mirror symmetry

Problem statement with the stream function

We solve the problem with the object in the formulation( the same formalism as in the potential flow) using
the stream function (script p.117). Show that the 3 equations (1) can be written in the following form using
, the vorticity and the pressure p :

2 = (2)

2 = (3)
y x x y
 2 2
1 2 2 2
p = 2 2
, (4)
2 x y xy

where the vorticity and the stream function are given by,

vx = vy = = ( v)z (5)
y x
In the problem we are only interested in the scalar functions and , the uncoupled equation (4) need not
to be considered. Additionally because of the mirror symmetry of the problem, the functions and satisfy
(x, y) = (x, y) and (x, y) = (x, y), thus we only need to simulate the upper half of the system.

Asymptotic Solution- Stream function far from the rectangular object

Far from the object we can easily determine the velocity field using the translation symmetry in the x-
coordinate. The velocity field must be of the form v0 = vx (y)ex .Thus the continuity equation is satisfied
trivially. The y-component of the Navier-Stokes equation yields y = 0. Verify that the x-component of the
Navier-Stokes equation is solved for non-slip boundary conditions at y = H by

p H 2
  y 2 
v0 (y) = v0 1 ex with v0 = (6)
H l 2
where p/l is the pressure gradient. The respective stream function and vorticity yield,
 
1  y 2
0 (y) = v0 y 1 , 0 (y) = 2v0 y/H 2 (7)
3 H

Discretization of the problem

For numerical purposes, the problem must be discretized. For simplicity we choose a regular discretization
with discretization length h, see Fig. 2. Additionally it is practical to work with dimensionless quantities, so we
define = , = , y = y/h,
v0 h v0 /h x = x/h and v = v/v0 , then the Eq. (2-3) are rescaled in the following way:

2 =
" !#

2 R
= 0 (9)
y x
x y

where R = v0h is the Reynolds number for our discretization. Note that the differential operators are also
affected with the new spatial coordinates x
and y. We also use the Finite Differences method to numerically
solve the above equations; We use symmetrical form for the derivatives,

2 f = fi(j+1) + fi(j1) + f(i+1)j + f(i1)j 4fij 
f = f(i+1)j f(i1)j /2 (10)
ij x
If one writes the discretized Eq. (8-9) in matrix form, one obtains in multi-index notation,

A(ij)(i0 j 0 ) (i0 j 0 ) = q(ij) and B(ij)(i0 j 0 ) (i0 j 0 ) = r(ij) (11)

The matrix A is constant and the vector q is given by the vorticity and the boundary conditions at .

The matrix B depends on and the vector r is equal to 0, except for a few components which determine the
boundary conditions. The solution is obtained in an iterative process:

and (for
1. Initialize 0 and0 );
example, by
= q and obtain the new stream function ;
2. Solve A

3. Solve B = r and get new vorticity ;

4. Repeat points 2 and 3 until (approximate) convergence

Figure 2: Discretization of the problem, only the upper half is shown; W = wh, T = th.

Boundary conditions
The stream function has been chosen to have the value 0 on the axis of symmetry y = 0. The edges of the
object, together with the symmetry axis, form a single stream line, hence on the object edges must be equal
to 0. The upper edge at y = H is also a stream line, since we can also use the Dirichlet boundary condition
(y = H) = 0 (y = H) = (2v0 H)/3. For the left and right edges, we can write the Neumann boundary
condition x = 0 (translation invariance).

The vorticity must vanish on the axis of symmetry on the basis of the symmetry, ie. (y = 0) = 0. For the

left and right edges, we use Neumann boundary conditions x = 0 (translation invariance). Slightly less obvious
are the conditions for the stream boundaries at the object and at the upper edge. For this reason the non-slip
condition for viscous liquid with Taylor expansion of the stream function comes to the aid - for the upper edge
of the object , for example, the vorticity is

2 2
iw = = (12)
y2 iw


the expansion of the stream function gives

1 2
i(w+1) = iw + + + .... (13)
y iw 2 y2 iw

Let us now consider that the non-slip condition at the upper object edge states that y = vx = 0 , we get
 
iw = 2 i(w+1) iw (14)
Similar conditions are to be applied to other edges. For the corners of the object, a symmetrical form is
offered by the condition for the upper and the respective side edge, for the left corner
1 Note: a Dirichlet boundary condition using the asymptotic solution 0 could also be used.

aw = a(w+1) + (a1)w 2aw
inside the object - a
We need additional (fictitious) boundary conditions for the apparent variable ,
Dirichlet condition that fixes the value to 0 satisfies this.

The skeleton of the sample program

The program consists of 4 modules written in the programming language Python, contains the iterative solution for the coupled equations; here are determined all physical parameters of the problem ; contains the skeletons of the functions that implement the eq. (11), which must be
supplemented in this problem; includes everything for the graphical representation of the solution.

The Python libraries NumPy, SciPy and Matplotlib are used for this purpose. To start the program, down-
load the 4 files from the continuum mechanics web page, save them to a directory, and execute the command
python in the directory.

Complete the missing lines with discretized equations and the boundary conditions for the upper, lower,
left, and right edges, and test the program in debug_mode (i.e., without the object). 2
Complete the boundary conditions for the object, switch off debug_mode, and make test simulations for
the problem with H = 4W , W = T with (physical) Reynolds number Re 2W v0 / equal to 3 and 24.

Increase the precision of the self-consistent calculation, given by variable aimed_precision in the module (to about 108 ).
Investigate the effect of the boundary conditions - increase the distance in front of the object upstream
(ah F ) together with the length to the right edge (L F T ) and find suitable F and L.

Investigate the effect of discretization by increasing the node points density (points_per_unit_length
in the module ).

2 Note: the analytical solution in Eq.(7) for problem without object is not an exact solution of the discretized problem.