Biringen S., Chow C.-Y. An Introduction to Computational Fluid Mechanics by Example

Подождите немного. Документ загружается.

172 VISCOUS FLUID FLOWS

0 0.2 0.4 0.6 0.8 1

0.2

0.4

0.6

0.8

Velocity (m/s)

Height (m)

FIGURE 3.5.2 Velocity distribution at different times (in seconds). +, T = 25,000;

◦

T = 50,000; , T = 75,000; ♦, T = 100,000; ∇, T = 125,000.

Problem 3.10 Find the velocity distribution in the channel described in Pro-

gram 3.4 with the upper plate replaced by one oscillating at the speed u

sin ωt,

where u

= 1m/sandω = 1/1000 s

−1

Problem 3.11 In approximating the differential equation (3.5.1) by the ﬁnite-

difference equation (3.5.3), we used forward difference in time and central differ-

ence in space, so that the truncation error of the approximation is O(τ , h

).An

improved explicit scheme with truncation error of O(τ

, h

) can be constructed

by also using central difference in time and by replacing the second u

i, j

on the

right side of (3.5.3) by the time-average (u

i, j −1

i, j +1

)/2. Thus,

i, j +1

−u

i, j −1

2τ

= ν

i−1, j

−u

i, j −1

−u

i, j +1

i+1, j

(1 + 2R)u

i, j +1

= 2R(u

i−1, j

i+1, j

) +(1 −2R)u

i, j −1

(3.5.8)

This is called the DuFort–Frankel formula, which involves three time levels, as

does the formula (2.10.13) derived for hyperbolic differential equations.

Show that (3.5.8) is an unconditionally stable numerical scheme.

IMPLICIT METHODS FOR SOLVING PARABOLIC PARTIAL DIFFERENTIAL EQUATIONS 173

Let us take a close look at the explicit formula (3.5.3). The solution at the

grid point (i, j +1) is computed, using this formula, from the solutions evaluated

at three grid points (i−1, j), (i, j ),and(i + 1, j). These values, in turn, are

computed from solutions in their neighborhood at the previous time step. In this

way we can trace out the region of dependence of the point (i, j + 1),whichis

conﬁned between the two dashed lines shown in Fig. 3.5.1. This means that the

disturbance created at any other height in the ﬂuid reaches the height y

with a

ﬁnite speed h/τ . This contradicts the real situation in an incompressible ﬂuid, in

which a disturbance at any point is felt immediately by all parts throughout the

ﬂuid. Thus, to improve the accuracy of (3.5.3), we may reduce the size of τ or

the value of R. In so doing the dashed lines will approach the horizontal grid line

passing through (i, j + 1) and, in the meantime, more time steps will be needed

in the computation to reach the same time level. The improved accuracy in the

explicit method is therefore obtained at the expense of an increased amount of

computer time.

3.6 IMPLICIT METHODS FOR SOLVING PARABOLIC PARTIAL

DIFFERENTIAL EQUATIONS—STARTING FLOW IN A CHANNEL

The deﬁciency associated with the explicit methods, that the solution computed

at one point is not affected immediately by the conditions at all other points in the

ﬂuid, can be avoided by devising an alternative numerical scheme for solving the

same diffusion equation (3.5.1). If we still use centered difference in space but

use backward, instead of forward, difference in time, a ﬁnite-difference equation

is obtained at (i, j) of the form

i, j

−u

i, j −1

) =

i−1,j

−2u

i, j

i+1,j

) (3.6.1)

It becomes, after regrouping,

i−1, j

−(1 + 2R)u

i, j

+Ru

i+1, j

=−u

i, j −1

(3.6.2)

where R = ντ/h

is the same dimensionless parameter as that deﬁned in

Section 3.5.

The slight modiﬁcation in approximating the time derivative causes a radical

change in the procedure for obtaining a solution. Suppose the solution at t = t

is to be computed based on the solution known at the previous time step t =

j−1

; (3.6.2) shows that every three neighboring unknown values are interrelated

through this linear algebraic equation. Applying (3.6.2) at all grid points interior

to the boundaries at one time level gives a system of simultaneous equations that

can be solved for all the unknowns at that time instant. In this way the velocities

at different heights are not independent of one another; a change at one point

will be felt immediately by all other points. Thus, this numerical scheme is more

sound than the explicit scheme on physical grounds. By using the numerical

174 VISCOUS FLUID FLOWS

scheme the solution can no longer be computed explicitly as before, so (3.6.2)

is called a formula for the implicit method.

The computational stability of (3.6.2) can be examined again with von Neu-

mann’s stability analysis by assuming the form already shown in (3.5.4) for the

numerical solution. It can easily be veriﬁed that the resulting relationship from

that analysis is

1 + 2R(1 −cos kh)

j−1

= λU

j−1

(3.6.3)

As cos kh varies from −1to+1, the value of the ampliﬁcation factor λ changes

from 1/(1 +4R) to 1 and can never exceed 1. Therefore, this numerical scheme

is stable for all positive values of R.

Although for computational stability there is no restriction on the magnitude

of R as long as it is positive, a smaller value of R results in a more accurate

numerical solution. The reason for this is that after multiplying (3.6.1) through by

τ , the truncated higher-order terms on the right-hand side are all multiplied by R.

We now apply the implicit method to solve a problem concerning the develop-

ment of a channel ﬂow caused by the application of a constant pressure gradient.

The initially stationary incompressible ﬂuid contained between two parallel inﬁ-

nite plates is set in motion by a suddenly imposed pressure gradient dp/dx along

the channel. Simpliﬁed for the present geometry, the equation of motion (3.1.7)

becomes

∂u

∂t

=−

+ν

∂

∂y

(3.6.4)

If the distance between plates is 2L and the origin of the coordinate system is

placed at the middle of the channel, the boundary and initial conditions are

u = 0aty =±L for all t (3.6.5)

u = 0att = 0for−L ≤ y ≤ L (3.6.6)

We know that as time increases, the velocity proﬁle will approach its steady-

state parabolic distribution

=−

2μ

−y

) (3.6.7)

which is a particular solution to (3.6.4) satisfying the boundary conditions (3.6.5).

By introducing the dimensionless variables

T =

/ν

, Y =

, U = u



−

2μ



(3.6.8)

and a dimensionless velocity difference

W = (u

−u)



−

2μ



= (1 −Y

) − U (3.6.9)

IMPLICIT METHODS FOR SOLVING PARABOLIC PARTIAL DIFFERENTIAL EQUATIONS 175

the governing equation (3.6.4) is simpliﬁed to

∂W

∂T

∂

∂Y

(3.6.10)

with boundary and initial conditions

W = 0atY =±1forallT (3.6.11)

W = 1 −Y

at T = 0for− 1 ≤ Y ≤ 1 (3.6.12)

The implicit numerical scheme for solving (3.6.10) is, according to (3.6.2),

i−1, j

−(1 +2R)W

i, j

+RW

i+1, j

=−W

i, j −1

(3.6.13)

where R = τ/h

, τ and h being the interval sizes in the nondimensionalized time

and space coordinates, respectively. Equation (3.6.13) is, in fact, a special case

of the ﬁnite-difference equation (2.2.10) with constant coefﬁcients C

on the

left-hand side, so that, after it is applied to all interior grid points, the resulting

tridiagonal system of simultaneous equations can be solved by calling the sub-

routine

TRID constructed in Program 2.1. Comparing (3.6.13) with the standard

form (2.2.10) reveals that the coefﬁcients are

= R, C

=−1 −2R, C

= R, C

=−W

i, j −1

The subscripts j and j −1 in (3.6.13) can be dropped in the computer program

if we compute the coefﬁcients C

using the values of W

at hand, call TRID to

get the solution at the next time step that will be stored under the same name

, and then advance in time by repeating the same process. Moreover, because

the ﬂow is symmetric about the x axis, it sufﬁces to ﬁnd a solution only in the

upper half of the channel. Accordingly, the conditions (3.6.11) and (3.6.12) are

modiﬁed to

W = 0atY = 1forallT (3.6.14)

∂W

∂Y

= 0atY = 0forallT (3.6.15)

W = 1 −Y

at T = 0for0≤ Y ≤ 1 (3.6.16)

The derivative boundary condition (3.6.15) indicates that a ﬁctitious grid point

is needed below Y = 0. Thus, to use subroutine

TRID, which solves n simulta-

neous equations, we divide the range 0 ≤ Y ≤ 1inton segments of equal width

h, denote the grid points by Y

, Y

, ..., Y

n+1

, respectively, and ﬁnally add a ﬁc-

titious point Y

at Y =−h. Equation (3.6.15) then becomes W

= W

,andthe

boundary conditions are exactly analogous to those for the temperature proﬁle

on the ﬂat-plate thermometer in Section 3.3. Following the discussion in that

section, before

TRID is called, the coefﬁcient matrix must be modiﬁed using the

176 VISCOUS FLUID FLOWS

following assignment statements:

← (C

) (3.6.17)

← 0 (3.6.18)

← (C

−C

n+1

) (3.6.19)

← 0 (3.6.20)

The dimensionless velocity U at any time level is calculated from (3.6.9) once

W becomes known for that time.

Velocity proﬁles are printed (Table 3.A.5) and plotted selectively in Program

3.5, computed for n = 20 and R = 0.4. In the plot the nondimensionalized steady-

state velocity distribution U

is added to compare with the proﬁles at other time

instants.

One thousand time steps have been marched in the computation of Pro-

gram 3.5 without encountering any instability, as expected of the implicit method.

Figure 3.6.1 displays the development of the ﬂow from a stationary state to the

steady-state parabolic velocity proﬁle as time approaches inﬁnity. Right after the

pressure gradient is applied, the ﬂow acceleration is large, but the acceleration

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

0.2

0.4

0.6

0.8

Dimensionless velocity

Dimensionless height

FIGURE 3.6.1 Velocity proﬁles. ×, T = 0;

◦

, T = 0.2; , T = 0.4; ♦; T = 0.6; ∇, T = 0.8;

, steady state.

IMPLICIT METHODS FOR SOLVING PARABOLIC PARTIAL DIFFERENTIAL EQUATIONS 177

diminishes with increasing time. At T = 1 the velocity at the center of the chan-

nel is 91% the steady-state value. For water (ν = 1 × 10

−6

/s) in a channel

2 cm high, this dimensionless time corresponds to 100 s in real time. The slow

readjustment of the vorticity distribution is a characteristic of laminar ﬂows with

small kinematic viscosity.

The governing equation (3.6.10) and the boundary and initial conditions

(3.6.11) and (3.6.12) for W show that W (Y , T ) describes the deceleration of a

channel ﬂow with a parabolic velocity proﬁle when the pressure gradient used

to sustain the ﬂow is suddenly removed.

For the problem considered here, an analytical solution to (3.6.4) is possible

by assuming

u =−

2μ

−y

) +

∞



m=1

exp



−λ

νt



cos





The eigenvalues λ

are to be determined from the boundary conditions (3.6.5),

and the coefﬁcients A

are to be evaluated after substituting this expression into

the initial condition (3.6.6) and using the orthogonal relationships among cosine

functions.

The truncation error resulting from approximating the differential equation

(3.5.1) by the implicit formula (3.6.2) is O(τ, h

). Following a procedure very

similar to that used to obtain the DuFort-Frankel formula (3.5.8) for improve-

ment of the explicit method, we may derive an improved implicit scheme with

a truncation error of O(τ

, h

). Let us replace the differential equation by a

ﬁnite-difference equation at the midpoint (i, j −

) between grid points (i, j ) and

(i, j −1), as shown in Fig. 3.6.2. The time derivative in (3.5.1) is approximated

by the central-difference formula (2.2.8) to give



∂u

∂t



i, j −(1/2)

i, j

−u

i, j −1

2(τ/2)









(i + 1, j )

(i, j )

(i − 1, j )

(i − 1, j − 1)

(i, j − 1) (i + 1, j − 1)

(i, j + )

τ/2

FIGURE 3.6.2 Crank-Nicolson method.

178 VISCOUS FLUID FLOWS

The second-order spatial derivative at (i, j −

) is approximated by the average of

those at (i, j) and (i, j − 1), each of which is expressed in the central-difference

form. Thus, we have



∂

∂y



i,j−(1/2)



i−1,j

−2u

i, j

i+1,j

i−1,j−1

−2u

i,j−1

i+1,j−1



+O(h

)

Substituting these two approximations into (3.5.1) yields, after rearranging and

dropping terms O(τ

) and O(h

i−1, j

−2(1 + R)u

i, j

+Ru

i+1, j

=−Ru

i−1, j −1

−2(1 −R)u

i, j −1

−Ru

i+1, j −1

(3.6.21)

The form of this more accurate formula differs from the form of (3.6.2) in that

it involves the evaluation of u at three grid points instead of at one on the

right-hand side. The corresponding increase in computational work is merely

in the evaluation of the coefﬁcients C

, which is apparently an insigniﬁcant

portion of the total computational effort. The numerical technique using (3.6.21)

is called the Crank-Nicolson method. It can easily be shown that this is also an

unconditionally stable scheme.

Problem 3.12 For the unsteady ﬂow in an inﬁnitely long tube of circular cross

section along the x axis, the governing equation is

∂u

∂t

=−

∂p

∂x

+ν



∂

∂r

∂u

∂r



(3.6.22)

Derive an implicit scheme for this differential equation in a form analogous

to (3.6.2), examine its stability, and then use it to ﬁnd the solution for the ﬂow

caused by the sudden application of a constant pressure gradient along a tube of

radius A. The analytical solution of this problem can be found in Section 4.3 of

Batchelor (1967), expressed in terms of a Fourier-Bessel series.

Project for Further Study: Find the solution for a ﬂow in a tube of radius A

caused by the application of an alternating pressure gradient

∂p

∂x

=−K cos ωt

where K is the constant amplitude and ω is the circular frequency of the oscilla-

tion. For this problem it is more convenient to introduce dimensionless variables

T = ωt, R =

, U =

/4μ

so that the governing equation (3.6.22) becomes, with  = ωρ A

/μ,



∂U

∂T

= 4cosT +

∂

∂R

∂U

∂R

NUMERICAL SOLUTION OF BIHARMONIC EQUATIONS—STOKES FLOWS 179

Using the numerical scheme constructed in the previous problem, ﬁnd U (R, T )

for  = 0.1, 7, and 35 and compare them with the analytical solution plotted in

Fig. 1 of Chow and Lai (1972).

3.7 NUMERICAL SOLUTION OF BIHARMONIC

EQUATIONS—STOKES FLOWS

In the previous three sections, examples with special geometries were chosen so

that the nonlinear terms contained in the substantial derivative on the left side

of (3.1.7) disappeared. The same result can be achieved by requiring that the

Reynolds number of the ﬂow be small.

Assuming that the incompressible ﬂow is characterized by a length L and a

speed U , we may introduce the following dimensionless variables:



L/U

, (x



, y



, z



) =





, V



, p



ρU

In terms of the new variables, the Navier-Stokes equation (3.1.7) becomes



=−∇



∇

2

V (3.7.1)

where the primed operators are those in dimensionless form and Re = ρUL/μ is

the Reynolds number.

As stated earlier in Section 2.1, the Reynolds number (Re) measures the rel-

ative importance of inertial force and viscous force in a ﬂow. If Re is large, the

viscous force terms, or the second group on the right side of (3.7.1), become

small in comparison with the inertial force terms on the left. In this case the

viscous force is important only in a thin region on the surface of a body in

motion relative to the ﬂuid. This thin region is the boundary layer considered in

Sections 3.2 and 3.3. On the other hand, if Re is much smaller than unity, viscous

force has a dominating inﬂuence on the ﬂuid motion. By dropping the inertial

force term and returning to dimensional notation, we obtain for the steady state

the linear equation

μ∇

V = ∇p (3.7.2)

Flows for which Re  1 are called the Stokes ﬂows or creeping ﬂows.Theﬂow

of molasses and the ﬂow induced by the motion of a microorganism are examples.

Taking the curl of (3.7.2) gives

∇

ζ = 0 (3.7.3)

in which ζ = ∇ ×V is the vorticity vector. Alternatively, divergence of (3.7.2)

yields, remembering that ∇ · V = 0 for incompressible ﬂuids,

∇

p = 0 (3.7.4)

180 VISCOUS FLUID FLOWS

This states that the pressure in an incompressible Stokes ﬂow satisﬁes the Laplace

equation.

Here we consider the special case of two-dimensional Stokes ﬂows parallel to

the x-y plane. If the stream function introduced in Section 2.1 is used so that

u =

∂ψ

∂y

and v =−

∂ψ

∂x

(3.7.5)

the only nonvanishing vorticity component ζ in the z direction can be written,

according to (2.8.16), as

∇

ψ =−ζ (3.7.6)

and the vector equation (3.7.3) reduces to a scalar equation

∇

ζ = 0 (3.7.7)

Combining (3.7.6) and (3.7.7) results in

∇

ψ = 0 (3.7.8)

where the operator is deﬁned as

∇

=∇

∇

∂

∂x

∂

∂x

∂y

∂

∂y

(3.7.9)

Because the Laplace equation is also called the harmonic equation, (3.7.8) is

usually referred to as the biharmonic equation.

Numerical solution of the biharmonic equation (3.7.8) seems to be simple. It

may be treated as two simultaneous second-order equations (3.7.6) and (3.7.7) to

each of which the numerical methods developed in Section 2.8 for solving elliptic

equations may be applied. However, difﬁculties arise in specifying the boundary

conditions for the vorticity equation (3.7.7). On the boundaries of a ﬂow region

we can only prescribe the velocity components and/or some of their derivatives

or, equivalently, the stream function and/or some of its derivatives. Vorticity is to

be computed from the velocity ﬁeld, but it cannot be speciﬁed at the boundaries

before the problem is solved. Thus, the numerical solution of the biharmonic

equation is not as simple as it ﬁrst appears to be, and special techniques are

needed to utilize the numerical methods, devised originally for solving Poisson

and Laplace equations, to solve the present boundary-value problem.

An example with an especially simple geometry is used here to demonstrate

how such a problem is handled. A square cavity DABC is shown in Fig. 3.7.1

within which a steady ﬂuid motion is generated by sliding an inﬁnitely long plate

lying on top of the cavity. Suppose that all variables are normalized so that the

size of the cavity is 1 ×1 and the sliding velocity is 1 in the positive x direction,

and that the Reynolds number is so low that the induced motion can be classiﬁed

as a Stokes ﬂow.

NUMERICAL SOLUTION OF BIHARMONIC EQUATIONS—STOKES FLOWS 181

(i, m − 2)

j = m

j = 1

i = m

i = 1

(i +1, m −1)

(i, m)

(i −1, m −1)

(i, m − 1)

ψ = 0

ψ = 0,

∂ψ

∂y

= 0

ψ = 0,

∂ψ

∂y

= 1

ψ = 0

∂ψ

∂x

= 0

∂ψ

∂x

= 0

FIGURE 3.7.1 Cavity ﬂow caused by a moving plate.

If no ﬂuid is squeezed out of the cavity below the moving plate, the ﬂuid

motion forms closed paths within the cavity. The surfaces DA, AB, BC ,and

CD are then segments of the bounding streamline designated by ψ = 0. This is

equivalent to specifying that the velocities normal to these four surfaces are all

zero. To require that the tangential velocity be vanishing on all surfaces except

on the top plate where it is 1, we obtain four additional boundary conditions

in terms of derivatives of ψ, as indicated in Fig. 3.7.1, along the four sides. In

summary, the boundary conditions are

ψ = 0and

∂ψ

∂x

= 0onDA (3.7.10)

ψ = 0and

∂ψ

∂y

= 0onAB (3.7.11)

ψ = 0and

∂ψ

∂x

= 0onBC (3.7.12)

ψ = 0and

∂ψ

∂y

= 1onCD (3.7.13)

This problem is similar to the one considered in Section 7.2 of Greenspan (1974),

except in that analysis the inertial terms are retained. Furthermore, we will use

a different procedure to compute the vorticity at the boundaries.