Fluid Mechanics 4e Solutions - Kundu Cohen

Chapter 1 Page 1 of 7

11.1 Show that the stability condition for the explicit scheme (11.10) is the condition (11.26).
Solution
The stability of the explicit scheme
T
i
n+1
= T
i
n
T
i +1
n
T
i1
n
( )+ T
i +1
n
2T
i
n
+ T
i1
n
( )

with = u
t
2x
, = D
t
x
2
, may be examined using the von Neumann method. The condition
for stability is discussed in the text, and is given by g
n+1
g
n
1, where

g
n+1
g
n
= +
( )
e
ikx
+ 1 2
( )
+
( )
e
ikx
= + ( ) cos isin ( ) + 1 2 ( ) + ( ) cos + i sin ( )
= 1 2 ( )+ 2 cos i2 sin
= 1 2 1 cos ( ) i2 sin
= 1 4sin
2
2
( )
i2 sin

for any value of the wavenumber k, where = kx. Therefore, the stability condition reduces to
1 4 sin
2
2 ( )
( )
2
+ 2sin ( )
2
1.

11.2 For the heat conduction equation
T
t
D
2
T
x
2
= 0, one of the discretized forms is
sT
j +1
n+1
+ 1+ 2s ( )T
j
n+1
sT
j 1
n+1
= T
j
n

wheres= D
t
x
2
. Show that this implicit algorithm is always stable.
Solution
For this implicit scheme, the evolution of error satisfies
s
j +1
n+1
+ 1+ 2s
( )
j
n+1
s
j1
n+1
=
j
n
(1)
Consider a component of the error in the Fourier space

j
n
= g
n
k ( )e
ikx
j
,
where k is the wavenumber in Fourier space, and g represents the function g at time t
n
= t
n
. The
component at the next time level has a similar form

j
n+1
= g
n+1
k
( )
e
ikx
j
.
Substituting the above two expressions into the error equation (1), we have,
g
n+1
se
ikx
j +1
+ 1+ 2s ( )e
ikx
j
se
ikx
j 1
[ ]
= g
n
e
ikx
j
,

or

g
n+1
g
n
=
1
se
ikx
+ 1+ 2s ( ) se
ikx
[ ]
=
1
1+ 2s 2scos
=
1
1+ 2s1 cos
( )
=
1
1+ 4ssin
2
2
( )

where = kx. Since 4ssin
2
2 ( ) is always positive, thus g
n+1
g
n
1, the scheme is
unconditionally stable.

11.3 An insulated rod initially has a temperature of T x,0 ( )= 0C 0 x 1 ( ). At t = 0 hot
reservoirs are brought into contact with the two ends, A T =100C ( ) x= 0 ( ) x and B ( 1 = )
T 0,t ( )= T ,t
:
. Numerically find the temperature T x 1,t ( )= 100C ( ) of any point in the rod. The
governing equation of the problem is the heat conduction equation
T
t
D
2
T
x
2
= 0. The exact
solution to this problem is
T
*
x
j
,t
n
( )=100
400
2m1 ( )
sin 2m1 ( )x
j
[ ]
expD 2m1 ( )
2
2
t
n
[ ]
m=1
NM

where NM is the number of terms used in the approximation.
(a). Try to solve the problem with the explicit FTCS (forward time, central space) scheme. Use
the parameter s= D
t
x
2
=0.5 and 0.6 to test the stability of the scheme.
(b). Solve the problem with a stable explicit or implicit scheme. Test the rate of convergence
numerically, using the error atx=0.5.
Solution
(a). The FTCS scheme can be written as
T
i
n+1
= T
i
n
+ sT
i +1
n
2T
i
n
+ T
i 1
n
( )

where s= D
t
x
2
. In the program, a uniform grid spacing and constant time step are used, and
their values are x=130 (with 31 grid points in the domain) and t =1500 (with 500 time
steps reaching the final time of 1 second), respectively.
The exact solution is evaluated with 10 terms in the summation. It is found that 5 terms
would give almost the same solution.
The case with s=0.5 is first computed and the results are shown in figure 1. The temperature
distributions at four different times are plotted. It is shown that the solution is stable, and agrees
with the exact solution.

The results for the case with s=0.6 are presented in figures 2. It is observed that the solution
is not stable. Oscillations with ever larger amplitudes occur. Eventually the solution blows up.
The stability condition for this explicit scheme is s 0.5.

(b) Using the FTCS scheme, the rate of the convergence is numerically evaluated. In the
evaluation, the time step is fixed at
t 10 =
7
, and the grid spacing is reduced from

D= 1, the stability x= 0.1 to 0.00625. These setting are chosen such that, when the constant

condition, , is always satisfied. The differences between the numerical solution and the
exact solution (error) at
s 0.5
x= 0.5 and t = 0.1 are listed in the table below and are plotted in figure
3 as a function of grid spacing x. It is observed that the rate of the convergence (the slope of
the curve fitted line) is approximately 2.0.

Figure 3. Convergence of the numerical solution.

Table 1. Convergence test.
= T
num
T
exact
x Error

0.10000 1.2062E-02
0.050000 2.4300E-03
0.025000 5.7286E-04
0.012500 2.5716E-04
0.0062500 5.6250E-05

The FORTRAN source code for this program is listed below.
c****************************************************************
c
c pr ogr amsol ves t he unst eady heat conduct i on pr obl em
c usi ng an expl i ci t FTCS scheme
c
c****************************************************************
pr ogr amFTCS
i mpl i ci t none
i nt eger nx, nt , i , n, m
r eal s, dx, dt
par amet er ( nt =500, nx=30)
par amet er ( s=0. 5)

r eal *8 t ( 0: nt , 0: nx) , t s( 0: nt , 0: nx)
r eal *8 pi , t i me, x, D
c
open( uni t =12, f i l e=' out put ' )
c
pi = 3. 14159265358979d0
dt = 1. 0/ nt
dx = 1. 0/ nx
c
c i ni t i al condi t i ons
do i =1, nx- 1
t ( 0, i ) = 0
enddo
c
c boundar y condi t i ons
do n=0, nt
t ( n, 0) = 100
t ( n, nx) = 100

enddo
c
c numer i cal sol ut i on wi t h FTCS
do n=0, nt - 1
do i =1, nx- 1
t ( n+1, i ) = t ( n, i ) + s*( t ( n, i - 1) - 2. 0*t ( n, i ) +t ( n, i +1) )
enddo
enddo
c
c exact sol ut i on
D = s*dx**2/ dt
do n=1, nt
t i me = n*dt
do i =0, nx
x = i *dx
t s( n, i ) = 100. 0
do m=1, 10
t s( n, i ) = t s( n, i ) - ( 400. d0/ ( 2*m- 1) / pi ) *
& si n( ( 2*m- 1) *pi *x) *exp( - D*( ( 2*m- 1) *pi ) **2*t i me)
enddo
enddo
enddo

c
c pr i nt - out
do n=50, nt , 50
do i =0, nx
wr i t e( 12, ' ( 5e12. 5) ' )
& n*dt , i *dx, t ( n, i ) , t s( n, i ) , abs( t ( n, i ) - t s( n, i ) )
enddo
enddo
c
st op
end

11.4 Derive the weak form, Galerkin form, and the matrix form of the following strong problem:
Given functions D x , and constants g, h, find u x ( ), f x ( ) ( ) such that
[D(x)u
, x
]
,x
+ f(x) = 0 on = 0,1 ( )
,
with
u

and

0 ( )= g u
,x
(1) = h
.
Write a computer program solving this problem using piecewise linear shape functions. You may
set D g f x x = 1, = 1 h= 1 , and ( ) sin(2 = )
. Check your numerical result with the exact
solution.
Solution
Define the functional space and the variational space for the trial solutions,
S = u x ( ) uH
1
,u 0 ( ) = g
{ }
and V = w x ( ) wH
1
, w0 ( ) = 0
{ }
,
respectively. Multiply the governing equation by a function in the variational space (wV), and
integrate the product over the domain = 0,1 ( )
,


. [D(x)u
, x
]
, x
+ f ( )wdx
0
1
= 0
Integrating by parts, we have

Du
,x
w
, x ( )
dx
0
1
+ Du
, x
w
[ ]
0
1
+ f x
( )
wdx
0
1
= Du
,x
w
,x
( )
dx
0
1
D1 ( )hw1 ( )+ f x ( )wdx
0
1
= 0

where the boundary conditions, u
,x
(1) = h and w 0 ( ) = 0, are applied. Therefore, the weak form
can be stated as:
Given functions , and constant h, find u D x ( ), f x ( ) S such that for all wV,
. Du
, x
w
,x
( )dx
0
1
= f x ( )wdx hD1 ( )w1 ( )

0
1
We can construct finite-dimensional approximations to S and V, which are denoted by

S
h
and V
h h
, respectively. We can also define a new function g x ( )
h
( ) such that g 0 = g. The
Galerkin form of the problem can be expressed as,
Find u , where v , such that for allw
h
= v
h
+ g
h h
V
h h
V
h
,
a w
h
,v
h
( )
= hD 1 ( )w
h
1 ( ) + w
h
, f
( )
a w
h
,g
h
( )
,
where a w and . ,v ( ) = Dw
, x
v
,x
( dx
0
1
)
N
w, f ( ) = wf ( )dx
0
1
Next, we construct the finite-dimensional variational space V explicitly using the shape
functions,
h
A
x ( ), A=1,2,...,n,
. w
h
= c
A
N
A
x ( )
A=1
n
The function g can be expressed with an additional shape function

h
N
0
x ( ) ( ), N
0
0 ( ) = 1
g
h
x ( ) = gN
0
x ( ) , such that g
h
0 ( ) = g .
With these definitions, the approximate solution can be written as
. v
h
x ( )= d
A
N
A
x ( )
A=1
n
Substitution of the approximate solution into the Galerkin formulation yields

a c
A
N
A
A=1
n
, d
B
N
B
B=1
n

= hD1 ( ) c
A
N
A
1 ( )
A=1
n
+ c
A
N
A
A=1
n
, f

a c
A
N
A
A=1
n
, gN
0

.
Since the coefficients c are arbitrary, the above equation reduces to
A
for A=1,2,,n. d
B
a N
A
,N
B
( )
B=1
n
= hD1 ( )N
A
1 ( )+ N
A
, f ( ) ga N
A
,N
0
( )
Therefore, the matrix form of the problem can be written as
Kd= F ,

where
K = K
AB
[ ]
, F = F
A
{ }, d= d
B
{ }
, K
AB
= DN
B, x
N
A, x
( )dx
0
1
F
A
= hD1 ( )N
A
1 ( )+ N
A
f ( )dx
0
1
g DN
0, x
N
A, x
( )dx
0
1
.
Set D= 1, g= 1, and h= 1 f (x) = sin(2x)
,
the exact solution for this problem is found to
be
u x ( ) =
sin 2x ( )
4
2
1+
1
2

x+1.
In the numerical solution, we have the stiffness matrix,
, K
AB
= N
B,x
N
A,x
( )dx
0
1
More explicitly, with piece-wise linear shape functions, most of the elements in the stiffness
matrix are zero, except three non-zero ones (per equation),
K
AA
= N
A, x
2
dx
0
1
= N
A, x
2
dx+ N
A, x
2
dx
x
A
x
A+1
x
A1
x
A
=
1
h
2
dx+
1
h
2
dx
x
A
x
A+1
x
A1
x
A
=
2
h
,
K
A1, A
= K
A+1, A
= N
A, x
N
A1, x
dx
x
A1
x
A
=
1
h
1
h

dx
x
A1
x
A
=
1
h
,
and K
nn
= N
A, x
2
dx
x
n1
x
n
=
1
h
2
dx
x
A1
x
A
=
1
h
.
Therefore, the stiffness matrix is evaluated as,
K =
2 h 1h 0 ... 0 0
1h 2 h 1h ... 0 0
0 1h 2 h ... 0 0
. . . . . .
0 0 0 ... 2 h 1h
0 0 0 ... 1h 1h

.
It is a tridiagonal matrix.
In evaluating the force vector, we first interpolate the function by , thus f x ( ) = f
B
N
B
x ( )
B=0
n

F
A
= N
A
1
( )
+ f (x)N
A
dx
0
1
N
0,x
N
A,x
( )
dx
0
1
=
An
+ f
B
B=0
n
N
B
N
A
dx
0
1
+
1
h
A1
=
1
6
f
A1
+
2
3
f
A
+
1
6
f
A+1

h+
1
h
, A=1
1
6
f
A1
+
2
3
f
A
+
1
6
f
A+1

h, A= 2,...,n1
1+
1
6
f
A1
+
1
3
f
A

h, A= n

The program for the numerical solution is attached below. The numerical solution uses n= 50
uniform elements. The comparison of the numerical solution with the exact one is presented in
figure 4. The agreement is almost perfect.

-0.2
0
0.2
0.4
0.6
0.8
1
0 0.2 0.4 0.6 0.8 1
fem solution
exact solution
u
x

Figure 4. Comparison of FEM solution with the exact solution.
The FORTRAN program for this problem is attached below.
c
c sol vi ng a t r i di agonal syst emwi t h TDMA scheme
c a( i ) *v( i ) =c( i ) *v( i - 1) +b( i ) *v( i +1) +d( i )
c
pr ogr amFEM1D
i mpl i ci t none
i nt eger n
par amet er ( n=50)

r eal *8 a( n) , b( n) , c( n) , d( n) , v( n) , vexact ( n)

i nt eger i , j
r eal *8 h, pi , p( n) , q( n) , t

pi = 3. 141592654

h = 1. d0/ n
c
c l oad t he st i f f ness mat r i x
do i =1, n
a( i ) = 2. d0/ h
b( i ) = 1. d0/ h
c( i ) = 1. d0/ h
enddo
a( n) = 1. d0/ h
b( n) = 0. d0
c( 1) = 0. d0
c
c l oad t he f or ce vect or
t = 2. d0*pi *h
do i =1, n- 1
d( i ) = h/ 6. d0*( si n( t *( i - 1) ) +4. d0*si n( t *i ) +si n( t *( i +1) ) )
enddo
d( 1) = d( 1) +1. d0/ h
d( n) = - 1. d0+h/ 6. d0*( si n( t *( n- 1) ) +2. d0*si n( t *n) )
c
c TDMA sol ver
p( 1) = b( 1) / a( 1)
q( 1) = d( 1) / a( 1)
do i =2, n
t = a( i ) - c( i ) *p( i - 1)
p( i ) = b( i ) / t
q( i ) = ( d( i ) +c( i ) *q( i - 1) ) / t
enddo
v( n) = q( n)
do i =n- 1, 1, - 1
v( i ) = p( i ) *v( i +1) +q( i )
enddo
c
c exact sol ut i on
do i =1, n
t = i *h
vexact ( i ) = si n( 2. d0*pi *t ) / ( 4. d0*pi **2) - ( 1+1/ ( 2*pi ) ) *t +1. d0
enddo
c
c pr i nt r esul t s
do i =1, n
wr i t e( 12, ' ( 3( e12. 5, A) ) ' ) i *h, ' , ' , v( i ) , ' , ' , vexact ( i )
enddo
c
st op
end

11.5 Solve numerically the steady convective transport equation, u
T
x
= D
2
T
x
2
, for 0 x1,
with two boundary conditions T 0 ( ) = 0 and T 1 ( ) 1, where u and D are two constants, =
(a) using the central finite difference scheme in equation (11.91), and compare with the exact
solution;

(b) using the upwind scheme (11.93); and compare with the exact solution.
Solution
(a) When the central difference scheme is used,
0.5R
cell
T
j +1
T
j 1
( )
= T
j +1
2T
j
+ T
j 1
( ),
where the cell Peclet number R
cell
= ux D and the grid spacing x=1n. Or
1 0.5R
cell
( )T
j +1
2T
j
+ 1+ 0.5R
cell
( )T
j 1
= 0,
with T
0
= 0 and T
n
= 1.
The numerical solution of the problem is presented in the figure 5. In the solution the grid
space is selected as x= 0.05 n= 20 ( ). Three values of the Peclet number are used,
R
cell
=1.0, 2.0, 3.0, in the computation. The scheme generates oscillatory solutions when
R
cell
> 2.0, near the boundary. The FORTRAN code for this case is attached below.

Figure 5. Numerical solution with the central difference scheme.
(b) When the first order upwind scheme is used,
R
cell
T
j
T
j 1
( )
= T
j+1
2T
j
+ T
j 1
( ),
or
T
j +1
2+ R
cell
( )T
j
+ 1+ R
cell
( )T
j 1
= 0,
with T
0
= 0 and T
n
= 1.

The solution of the problem is presented in the figure 6. In the solution the grid space is selected
as . Two values of the Peclet number are used, x= 0.05 n= 20 ( ) R
cell
=1.0, 3.0. The scheme is
stable, and the solution does not generate oscillations in the boundary layer. However, the
solution is quite diffusive.

Figure 6. Numerical solution with the first order upwind scheme.

c
c sol vi ng a t r i di agonal syst emwi t h TDMA scheme
c a( i ) *v( i ) =c( i ) *v( i - 1) +b( i ) *v( i +1) +d( i )
c
pr ogr amhw5a
i mpl i ci t none
i nt eger n
par amet er ( n=20)

r eal *8 a( 0: n) , b( 0: n) , c( 0: n) , d( 0: n) , v( 0: n)

i nt eger i
r eal *8 h, p( 0: n) , q( 0: n) , t , r cel l

r cel l = 3. d0
h=1. d0/ n
c
c st i f f ness mat r i x
do i =1, n
a( i ) = - 2. d0
b( i ) = - ( 1. d0- 0. 5*r cel l )
c( i ) = - ( 1. d0+0. 5*r cel l )

enddo
a( 0) = 1. d0
b( 0) = 0. d0
c( 0) = 0. d0
a( n) = 1. d0
b( n) = 0. d0
c( n) = 0. d0
c
c f or ce vect or
do i =0, n- 1
d( i ) = 0. d0
enddo
d( n) = 1. d0
c
c TDMA sol ver
p( 0) = b( 0) / a( 0)
q( 0) = d( 0) / a( 0)
do i =1, n
t = a( i ) - c( i ) *p( i - 1)
p( i ) = b( i ) / t
q( i ) = ( d( i ) +c( i ) *q( i - 1) ) / t
enddo
v( n) = q( n)
do i =n- 1, 0, - 1
v( i ) = p( i ) *v( i +1) +q( i )
enddo
c
c pr i nt r esul t s
do i =0, n
wr i t e( 12, ' ( 2( e12. 5, A) ) ' ) i *h, ' , ' , v( i )
wr i t e( *, ' ( 2( e12. 5, A) ) ' ) i *h, ' , ' , v( i )
enddo
c
st op
end

6. Code the explicit MacCormack scheme with the FF/BB arrangement for the driven cavity flow
problem as described in Section 5. Compute the flow field at Re=100 and 400, and explore
effects of Mach number and the stability condition (11.110).
Sample c-code:
/ /
/ / Aut hor s: Howar d H. Hu and Andy E. Per r i n
/ / Expl i ci t scheme usi ng FF/ BB MacCor mack met hod f or a dr i ven cavi t y pr obl em
/ /

#i ncl ude <st dl i b. h>
#i ncl ude <st r i ng. h>
#i ncl ude <st di o. h>
#i ncl ude <mat h. h>


#def i ne MI N( i , j ) ( ( i ) <( j ) ? ( i ) : ( j ) )

/ / gl obal s
i nt i t , i t 0, nt i me, npr i nt ;
doubl e cur t i me;
i nt i var =0;
i nt next =0;

i nt nx, ny, node;
doubl e Re, M, L, H, dx, dy, dt ;
doubl e *r oe_n, *r oeU_n, *r oeV_n, *r oe_s, *r oeU_s, *r oeV_s, *u, *v;
i nt nbd;
doubl e c1, c2, c3, c4, c5, c6, c7, c8, c9, c10, c11;
doubl e *var , umi n, umax;

voi d I ni t i al ( voi d)
{
i nt i , k;

nx = 101;
ny = 101;
Re = 400. 0;
M = 0. 05;
L = 1. 0;
H = 1. 0;
dx = L/ ( nx- 1) ;
dy = H/ ( ny- 1) ;
dt = 0. 8*M*dx/ sqr t ( 3. 0) ;
nt i me = ( i nt ) ( 80. 0/ dt ) ;
npr i nt = 25;

pr i nt f ( " dt =%13f nt i me=%i \ n" , dt , nt i me) ;

c1 = dt / dx;
c2 = dt / dy;
c3 = dt / ( dx*M*M) ;
c4 = dt / ( dy*M*M) ;
c5 = 4. 0*dt / ( 3. 0*Re*dx*dx) ;
c6 = dt / ( Re*dy*dy) ;
c7 = dt / ( Re*dx*dx) ;
c8 = 4. 0*dt / ( 3. 0*Re*dy*dy) ;
c9 = dt / ( 12. 0*Re*dx*dy) ;
c10 = ( - 2. 0) *( c5+c6) ;
c11 = ( - 2. 0) *( c7+c8) ;

node = nx*ny;
r oe_n = cal l oc( node, si zeof ( doubl e) ) ;
r oeU_n = cal l oc( node, si zeof ( doubl e) ) ;
r oeV_n = cal l oc( node, si zeof ( doubl e) ) ;
r oe_s = cal l oc( node, si zeof ( doubl e) ) ;
r oeU_s = cal l oc( node, si zeof ( doubl e) ) ;
r oeV_s = cal l oc( node, si zeof ( doubl e) ) ;
u = cal l oc( node, si zeof ( doubl e) ) ;
v = cal l oc( node, si zeof ( doubl e) ) ;
var = cal l oc( node, si zeof ( doubl e) ) ;


f or ( i =0; i <nx; i ++)
{
k = ny*i +( ny- 1) ;
u[ k] = 1. 0;
r oeU_n[ k] = u[ k] ;
r oeU_s[ k] = u[ k] ;
}

i t = 0;
i t 0 =0;
cur t i me = 0. 0;

}

voi d Next I t er ( voi d)
{
i nt i , j , k, k1, k2, k3, k4;

/ / node or der i ng
/ / k1=( i , j +1)
/ / k4=( i - 1, j ) k= ( i , j ) k3=( i +1, j )
/ / k2=( i , j - 1)

i t 0 = i t 0+1;
cur t i me = cur t i me + dt ;

/ / Pr edi ct or
f or ( i =1; i <nx- 1; i ++)
{
f or ( j =1; j <ny- 1; j ++)
{
k = i *ny+j ;
k1 = k+1;
k2 = k- 1;
k3 = k+ny;
k4 = k- ny;

r oe_s[ k] = r oe_n[ k] - c1*( r oeU_n[ k3] - r oeU_n[ k] )
- c2*( r oeV_n[ k1] - r oeV_n[ k] ) ;

r oeU_s[ k] = r oeU_n[ k]
- c1*( r oeU_n[ k3] *u[ k3] - r oeU_n[ k] *u[ k] )
- c2*( r oeV_n[ k1] *u[ k1] - r oeV_n[ k] *u[ k] )
- c3*( r oe_n[ k3] - r oe_n[ k] )
+ c10*u[ k] + c5*( u[ k3] +u[ k4] ) + c6*( u[ k1] +u[ k2] )
+ c9*( v[ k3+1] +v[ k4- 1] - v[ k3- 1] - v[ k4+1] ) ;

r oeV_s[ k] = r oeV_n[ k]
- c1*( r oeU_n[ k3] *v[ k3] - r oeU_n[ k] *v[ k] )
- c2*( r oeV_n[ k1] *v[ k1] - r oeV_n[ k] *v[ k] )
- c4*( r oe_n[ k1] - r oe_n[ k] )
+ c11*v[ k] + c7*( v[ k3] +v[ k4] ) + c8*( v[ k1] +v[ k2] )
+ c9*( u[ k3+1] +u[ k4- 1] - u[ k4+1] - u[ k3- 1] ) ;

}
}


/ / BC
f or ( j =1; j <ny- 1; j ++)
{
k=j ;
r oe_s[ k] = r oe_n[ k]
- 0. 5*c1*( - r oeU_n[ k+2*ny] +4. 0*r oeU_n[ k+ny] - 3. 0*r oeU_n[ k] )
- 0. 5*c2*( r oeV_n[ k+1] - r oeV_n[ k- 1] ) ;
k = ( nx- 1) *ny+j ;
- 0. 5*( - c1) *( - r oeU_n[ k- 2*ny] +4. 0*r oeU_n[ k- ny] - 3. 0*r oeU_n[ k] )
- 0. 5*c2*( r oeV_n[ k+1] - r oeV_n[ k- 1] ) ;
}
f or ( i =0; i <nx; i ++)
{
k = i *ny;
- 0. 5*c1*( r oeU_n[ k+ny] - r oeU_n[ k- ny] )
- 0. 5*c2*( - r oeV_n[ k+2] +4. 0*r oeV_n[ k+1] - 3. 0*r oeV_n[ k] ) ;
k = i *ny+( ny- 1) ;
- 0. 5*c1*( r oeU_n[ k+ny] - r oeU_n[ k- ny] )
- 0. 5*( - c2) *( - r oeV_n[ k- 2] +4. 0*r oeV_n[ k- 1] - 3. 0*r oeV_n[ k] ) ;
}

f or ( k=0; k<node; k++)
{
u[ k] = r oeU_s[ k] / ( 1. 0+r oe_s[ k] ) ;
v[ k] = r oeV_s[ k] / ( 1. 0+r oe_s[ k] ) ;
}

/ / Cor r ect or
f or ( i =1; i <nx- 1; i ++)
{
f or ( j =1; j <ny- 1; j ++)
{
k = i *ny+j ;
k1 = k+1;
k2 = k- 1;
k3 = k+ny;
k4 = k- ny;

r oe_n[ k] = 0. 5*( ( r oe_n[ k] +r oe_s[ k] ) - c1*( r oeU_s[ k] - r oeU_s[ k4] )
- c2*( r oeV_s[ k] - r oeV_s[ k2] ) ) ;

r oeU_n[ k] = 0. 5*( ( r oeU_n[ k] +r oeU_s[ k] )
- c1*( r oeU_s[ k] *u[ k] - r oeU_s[ k4] *u[ k4] )
- c2*( r oeV_s[ k] *u[ k] - r oeV_s[ k2] *u[ k2] )
- c3*( r oe_s[ k] - r oe_s[ k4] )
+ c10*u[ k] + c5*( u[ k3] +u[ k4] ) + c6*( u[ k1] +u[ k2] )
+ c9*( v[ k3+1] +v[ k4- 1] - v[ k3- 1] - v[ k4+1] ) ) ;

r oeV_n[ k] = 0. 5*( ( r oeV_n[ k] + r oeV_s[ k] )
- c1*( r oeU_s[ k] *v[ k] - r oeU_s[ k4] *v[ k4] )
- c2*( r oeV_s[ k] *v[ k] - r oeV_s[ k2] *v[ k2] )
- c4*( r oe_s[ k] - r oe_s[ k2] )
+ c11*v[ k] + c7*( v[ k3] +v[ k4] ) + c8*( v[ k1] +v[ k2] )
+ c9*( u[ k3+1] +u[ k4- 1] - u[ k4+1] - u[ k3- 1] ) ) ;

}
} / / end f or i

/ / BC
f or ( j =1; j <ny- 1; j ++)
{
k = j ;
r oe_n[ k] = 0. 5*( r oe_n[ k] +r oe_s[ k]
- 0. 5*c1*( - r oeU_s[ k+2*ny] +4. 0*r oeU_s[ k+ny] - 3. 0*r oeU_s[ k] )
- 0. 5*c2*( r oeV_s[ k+1] - r oeV_s[ k- 1] ) ) ;
k = ( nx- 1) *ny+j ;
- 0. 5*( - c1) *( - r oeU_s[ k- 2*ny] +4. 0*r oeU_s[ k- ny] - 3. 0*r oeU_s[ k] )
- 0. 5*c2*( r oeV_s[ k+1] - r oeV_s[ k- 1] ) ) ;
}
f or ( i =0; i <nx; i ++)
{
k = i *ny;
- 0. 5*c1*( r oeU_s[ k+ny] - r oeU_s[ k- ny] )
- 0. 5*c2*( - r oeV_s[ k+2] +4. 0*r oeV_s[ k+1] - 3. 0*r oeV_s[ k] ) ) ;
k = i *ny+ny- 1;
- 0. 5*c1*( r oeU_s[ k+ny] - r oeU_s[ k- ny] )
- 0. 5*( - c2) *( - r oeV_s[ k- 2] +4. 0*r oeV_s[ k- 1] - 3. 0*r oeV_s[ k] ) ) ;
}

f or ( k=0; k<node; k++)
{
u[ k] = r oeU_n[ k] / ( 1. 0+r oe_n[ k] ) ;
v[ k] = r oeV_n[ k] / ( 1. 0+r oe_n[ k] ) ;
}

}

i nt mai n ( i nt ar gc, char *ar gv[ ] )
{
i nt i ;

I ni t i al ( ) ;
i var = - 1;
f or ( i t =0; i t <nt i me; i t =i t +npr i nt )
f or ( i =0; i <npr i nt ; i ++) Next I t er ( ) ;

f r ee( r oe_n) ;
f r ee( r oeU_n) ;
f r ee( r oeV_n) ;
f r ee( r oe_s) ;
f r ee( r oeU_s) ;
f r ee( r oeV_s) ;
f r ee( u) ;
f r ee( v) ;
f r ee( var ) ;
r et ur n 0;
}
Problem 2.
p = f(L, a, , , , U)
p, L, a, , , , U; n = 7, dimensional parameters.
Primary dimensions: M, L, t; r = 3, primary dimensions. Now,
[p] = ML
1
t
2
, [L] = L, [a] = L, [] = ML
3
, [] = ML
1
t
1
, [] = t
1
, [U] = Lt
1
.
Selecting repeating parameters: a, , U, m = r = 3 repeating parameters.
Then, n m = 4 dimensionless groups will result.
Setting up the dimensional equations, we have:
1
a
C
1
C
2
U
C
3
L
C
4
L
C
1
(ML
3
)
C
2
(Lt
1
)
C
3
L
C
4
Therefore,
L : C
1
3C
2
+ C
3
+ C
4
= 0
M : C
2
= 0
t : 3C
3
= 0
Therefore, C
2
= 0 = C
3
and C
1
= C
4
. Thus,
1

_
L
a
_
C
4
. (1)
Next,
2
a
C
5
C
6
U
C
7
C
8
L
C
5
(ML
3
)
C
6
(Lt
1
)
C
7
_
ML
1
t
1
_
C
8
Therefore,
L : C
5
3C
6
+ C
7
C
8
= 0
M : C
6
+ C
8
= 0
t : C
7
C
8
= 0
Therefore, C
6
= C
8
, C
7
= C
8
, and C
5
= C
8
. Thus,
2

_

aU
_
C
8
. (2)
2
Next,
3
a
C
9
C
10
U
C
11
C
8
L
C
9
(ML
3
)
C
10
(Lt
1
)
C
11
_
t
1
_
C
12
Therefore,
L : C
9
3C
10
+ C
11
= 0
M : C
10
= 0
t : C
11
C
12
= 0
Therefore, C
10
= 0, C
11
= C
12
, and C
9
= C
12
. Thus,
3

_
a
U
_
C
9
. (3)
Next,
4
a
C
13
C
14
U
C
15
p
1
L
C
13
(ML
3
)
C
14
(Lt
1
)
C
15
_
ML
1
t
2
_
1
Therefore,
L : C
13
3C
14
+ C
15
1 = 0
M : C
14
+ 1 = 0
t : C
15
2 = 0
Therefore, C
15
= 2, C
14
= 1, and C
13
= 0. Thus,
4

_
p
U
2
_
. (4)
From 1-4, we get
p
U
2
C
_
L
a
_
C
4
_

aU
_
C
8 _
a
U
_
C
9
. (5)
where, C, C
4
, C
8
, C
9
are constants.
Therefore, we may write
p
U
2
= C
1
_
L
a
_
C
2
(Re)
C
3
(St)
C
4
. (6)
where, C
1
, C
2
, C
3
, C
4
are constants.
3
Problem 3.
Morgan and Young invoke the following assumptions: (1). A collar-like stenosis may be
approximated by a smooth, rigid, axisymmetric constriction in a long straight tube, (2).
The eect of the stenosis geometry is dominant and any inuence of wall distensibility is
negligible, (3). Blood can be treated as a Newtonian uid at the ow rates encountered
in the large arteries where stenosis commonly occur, (4). Blood ow is laminar, and (5).
Steady ow assumption is acceptable although the arterial blood ow is pulsatile.
The stenosis is shown in the following gure: The dimensional variables are designated by
Figure 1: Axisymmetric stenosis.
. The axial coordinate and velocity are z and u.

U is the centerline velocity. The radial
coordinate is r and the corresponding velocity component is v. The dimensionless variables
are: r = r/R
0
, z = z/R
0
, R =

R/R
0
, u = u/

U
0
, v = v/

U
0
, U =

U/

U
0
, and p = p/

U
0
2
,
where,

U
0
is the average velocity in the unobstructed tube, and p is the pressure and is
the density.
The usual continuity, axial momentum, and radial momentum equations for axisymmetric
ow in a tube are considered. The dimensionless momentum equations involve the Reynolds
number, Re = 2R
0

U
0
/, where is the uid viscosity. By integrating the axial momentum
equation over the cross section of the tube and by invoking the no slip condition, u = v = 0 at
the wall of the vessel, the integral momentum equation is developed. Next, by multiplying
the axial momentum equation by ru, and integrating over the cross section, the integral
energy equation is developed. It is now noted that in these two equations, the terms due to
the viscous component of the normal stress in the axial direction (
2
u/z
2
) are negligible.
Next, an important observation is made in regard to the pressure gradient p/z. It is
noted that the pressure gradient p/z is independent of the radial coordinate r for ow in
a straight tube, whereas for ow in a constricted tube the pressure gradient will, in general,
vary over the cross section. However, if the integral energy equation is multiplied by R
2
the
resulting integral involving the pressure gradient is equal to the corresponding integral in
the integral momentum equation for two special cases: (i) p/z is independent of r; (ii) the
velocity u is independent of r. Since these two conditions are approached in a constriction,
i.e. in the gradually varying initial portion of the constriction the pressure gradient is nearly
constant and in the rapidly converging portion the velocity prole tends to become attened,
we may let,
_
R
0
r
p
z
dr R
2
_
R
0
ru
p
z
dr. (1)
4
This approximation would enable the elimination of the pressure gradient term in both the
integral momentum and integral energy equations. Then it is possible to combine the integral
momentum and energy equations into a single equation in terms of axial velocity:
1
2
R
2

z
_
R
0
ru
3
dr

z
_
R
0
ru
2
dr =
2
Re
_
R
2
_
R
0
r
_
u
r
_
2
dr + R
_
u
r
_
R
_
. (2)
The equation (2) is subject to (i) u = U at r = 0; (ii) u = 0 at r = R; (iii) u/r = 0 at
r = 0; (iv)
3
u/r
3
= 0 at r = 0; and, (v) the condition that the net ow through any cross
section must be the same for any incompressible uid which may be expressed by:
_
R
0
r u dr =
1
2
. (3)
In the above, the condition (iii) is derived from a consideration of the forces on a cylindrical
element having its axis along the tube centerline. If the pressure and the inertial forces
are to be nite as the radius of the element approaches zero, the viscous force, which is
proportional to u/r, must approach zero. The condition (iv) is developed by eliminating
the pressure between the axial momentum and radial momentum equations and considering
the resulting equation as r approaches zero.
Next, it is noted that at high Reynolds numbers, the prole must allow a thin region of
high shear near the wall in the converging section with a relatively at prole in the core.
To accommodate these requirements, Morgan and Young construct a polynomial t which
permits the shear near the wall to become large while maintaining a at core ow. This t
is given by,
u = U for 0
r
R
, and, u = a + b
_
r
R
_
2
+ c
_
r
R
_
4
for
r
R
1, (4)
where a, b, and c are unknown coecients and is the value of (r/R) at the juncture sepa-
rating the at and polynomial parts of the prole. The unknown coecients are determined
from the no slip condition along with two compatibility conditions u = U and u/r = 0
at (r/R) = . The constraint (v) enables expressing in terms of R and U. Thus the
polynomial t prole for u is entirely in terms of U, r, and R. The prole is now introduced
into the equation (2). The resulting rst order, non-linear ordinary dierential equation is
numerically solved by assuming that Poiseuille ow prevails far upstream of the stenosis.
The solution provides the desired velocity proles. These are plotted by Morgan and Young.
The wall shear stress is evaluated from

w
=
_
u
r
_
R
_
1 +
_
_
2
_
, (5)
where

R
is the slope d

R/d z of the wall. The results are included in the paper by Morgan
and Young.
5
Problem 4.
In a pure pressure-gravity ow, the eects of friction are negligible compared with the eects
due to changes of the external pressure and the elevation in the gravity eld. Lengthwise al-
terations of speed, pressure, area, etc., are brought about by changes in p
e
and z. Changes in
p
e
may be brought about by : (i) active muscle tone; (ii) elastic constrictions, or sphincters,
as where veins pass from the abdominal cavity to the thoracic cavity, and in the pulmonary
system; (iii) weights, (iv) pressurizing cus, and (v) clamps etc.,. Changes in z are important
because (i) in the vertical position of a human being, the hydrostatic pressure exceeds the
venous and arterial pressure levels;(ii) during aircraft maneuvers; and, (iiI) during certain
phases of space ight when eective g may be greatly increased.
Retaining only the terms in d(p
e
+gz) in the table of inuence coecients for the tube law,
P
n
1, and, n =
3
2
, (1)
the governing equations of a pure pressure-gravity ow are written as,
_
1 S
2
_
1
d
dx
=
2
3
3
2
d
dx
, (2)
and,
_
1 S
2
_
1
S
2
dS
2
dx
= +
1
3
3
2
d
dx
, (3)
where,
=
(p
e
+ gz)
K
p
. (4)
Dividing equation (3) by equation (2) and integrating, for a pure pressure-gravity ow we
have,
=
A
A
=
_
u
u
_
1
= S
4
, (5)
where,

denotes the value of the quantity at S = 1. For the tube law given by equation (1),
the wave speed is known to be,
c =
_
nK
p
_1
2
, n =
3
2
. (6)
From equations(5) and (6),
c
c
=
_

3
4
= S
3
. (7)
Also, from Bernoullis theorem,
(p p
)
1
2
c
2
= 1 S
8
+
g (z
z)
c
2
. (8)
With equation (8), and by the elimination of from equation (3) using equation (5), Shapiro
develops,
1
3
3
2
d = S
4
(1 S
2
)dS
2
. (9)
6
Integration of equation (9) between limits and
corresponding to S and 1, gives

1
3
3
2
(
) =
1
3
_
S
6
1
_
1
4
_
S
8
1
_
. (10)
Shapiro provides graphs of equation(10). The graphs show that increasing values of drive
S towards unity and decreasing values of drive S away from unity. From these, it is
concluded that: (i) Choking occurs when the value of continually increases with either
S < 1 or S > 1; (ii) Continuous transition through S = 1 may be achieved by means of an
increase in until S = 1 is reached, with d/dx becoming zero exactly at S = 1 followed
then by a decrease of .
7
Problem 5.
To determine the volume ux for the ow with the Power-law model, consider an annu-
lar volume element of length L and thickness dr in the ow, as shown in the gure below:
r
dr
a
L
Figure 2: Representative volume element.
Let there be a constant pressure gradient, p/L, across the element of length L.
Force due to pressure gradient on the annular element = +
p
L
(2rL) dr. (1)
This must be balanced by shear force which is given by,
Shear force = (2dr) L
d
dr
(r). (2)
Therefore,
d
dr
(r) =
p
L
r (3)
=
p
L
r
2
+
c
2
(4)
where c is a constant. The stress is nite at axis r = 0. Therefore c = 0 and hence, from 4,
=
p
L
r
2
. (5)
For a power-law uid,
=
n
. (6)
Therefore =
_
p
L
r
2
_1
n
. (7)
In this problem, since velocity is a function of r only,
=
du
dr
. (8)
8
where u(r) is velocity component in r direction.
From equations 7 and 8,
du
dr
=
_
p
L
r
2
_1
n
. (9)
The no-slip condition at the wall requires that
u = 0 at r = a. (10)
By integrating equation 9 and using equation 10, one obtains
u =
_
p
2L
_
1/n
n
n + 1
_
a
n+1
n
r
n+1
n
_
. (11)
Now, the ux for ow, Q, is given by
Q =
_
R
0
2rudr. (12)
With integration by parts and invoking the no-slip conditions at r = a,
Q =
_
a
0
r
2
_
du
dr
_
dr. (13)
With equation 11,
Q =
_
a
0
r
2
_
p
L
r
2
_
1/n
dr (14)
=
_
p
2L
_
1/n
n
3n + 1
a
3n+1
n
(15)
When n = 1,
Q =
_
p
2L
_

4
a
4
(16)
=
_
a
4
8
_
p
L
a
3n+1
n
which agrees with equations (17.17) and is the Poiseuille formula.
9
Problem 6.
To determine the volume ux for the ow with Herschel-Bulkley model, we rst note that
=
n
+
0
,
0
and = 0 , <
0
(1)
There is yield stress and as a consequence the ow region includes plug ow in the core as
shown in gure Let the radius of the plug ow region be r
p
.
Core region
a
r
p
Velocity profile
Figure 3: Representative volume element.
For r r
p
, (r) =
0
= constant (2)
Therefore in the core, = 0. (3)
Therefore in the core,
du
dr
= 0. (4)
Therefore,
u = constant (5)
= u
p
(say) (6)
(7)
For a constant pressure gradient, (p/L), across an element of length L in the core of
region, from a force balance,
2r
p
L
0
=
p
L
r
p
2
L (8)
Thus, r
p
=
2
0
(p/L)
(9)
0
=
p
L
r
p
2
(10)
Outside the core region, with u = u(r), just as in problem 5,
du
dr
=
_
p
2L
_1
n
(r r
p
)
1
n
. (11)
10
By integration,
u =
_
p
2L
_
1/n
n
n + 1
(r r
p
)
n+1
n
+ C. (12)
where C is a constant of integration.
Now, at r = a, u = 0 (no-slip condition). Therefore,
u =
_
p
2L
_
1/n
n
n + 1
_
(a r
p
)
n+1
n
(r r
p
)
n+1
n
_
. (13)
We can also calculate the plug velocity by setting r = r
p
in equation 12. Thus,
u
p
=
_
p
2L
_
1/n
n
n + 1
(a r
p
)
n+1
n
. (14)
The volume ux may be calculated from,
Q = r
p
2
u
p
+
_
a
rp
2rudr. (15)
With equations 13 and 14, we can evaluate equation 15.
Let = r
p
/a. (16)
Then, after integration and considerable algebra,
Q =
_
p
2L
_
1/n
n
n + 1
a
3n+1
n
_
2
(1 )
n+1
n
+ (1 + )(1 )
2n+1
n
2n
3n + 1
(1 )
3n+1
n
2n
2n + 1
(1 )
2n+1
n
_
(17)
When
0
= 0, the Herschel-Bulkley model reduces to the Power-law model. Therefore, with
0
= 0 and = 0, equation 17 reduces to the result of problem 5.
11

Fluid Mechanics 4e Solutions - Kundu Cohen

Uploaded by

Copyright:

Available Formats

Fluid Mechanics 4e Solutions - Kundu Cohen

Uploaded by

Document Information

Original Description:

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Fluid Mechanics 4e Solutions - Kundu Cohen

Uploaded by

Copyright:

Available Formats

Chapter 1 Page 1 of 7

, and the grid spacing is reduced from

= f x ( )wdx hD1 ( )w1 ( )

We can construct finite-dimensional approximations to S and V, which are denoted by

The function g can be expressed with an additional shape function

Substitution of the approximate solution into the Galerkin formulation yields

corresponding to S and 1, gives

You might also like