Biringen S., Chow C.-Y. An Introduction to Computational Fluid Mechanics by Example
Подождите немного. Документ загружается.
2 FLOW TOPICS GOVERNED BY ORDINARY DIFFERENTIAL EQUATIONS
and the prescribed initial condition. Such a problem is called an initial-value
problem.
To solve the problem numerically, the axis of the independent variable is
usually divided into evenly spaced small intervals of width h whose end points
are situated at
t
i
= t
0
+ih, i = 0, 1, 2, ...
(1.1.2)
The solution evaluated at the point t
i
is denoted by x
i
. Thus, by using a numerical
method, the continuous function x(t) is approximated by a set of discrete values
x
i
, i = 0, 1, 2, ..., as sketched in Fig. 1.1.1. Since h is small and f is an analytic
function, the solution at any point can be obtained by means of a Taylor’s series
expansion about the previous point:
x
i+1
≡ x(t
i+1
) = x(t
i
+h)
= x
i
+h
dx
dt
i
+
h
2
2!
d
2
x
dt
2
i
+
h
3
3!
d
3
x
dt
3
i
+··· (1.1.3)
= x
i
+hf
i
+
h
2
2!
f
i
+
h
3
3!
f
i
+···
where f
n
i
denotes d
n
f
dt
n
evaluated at (x
i
, t
i
). f is generally a function of both
x and t, so that the first-order derivative is obtained according to the formula
df
dt
=
∂f
∂t
+
dx
dt
∂f
∂x
Higher-order derivatives are obtained by using the same chain rule.
t
0
t
1
t
2
x(t)
t
x
i+1
x
i
x
i−1
x
2
x
1
x
0
x
0
t
i
t
i+1
t
i−1
h
FIGURE 1.1.1 Numerical solution of an ordinary initial-value problem.
NUMERICAL SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS: INITIAL-VALUE PROBLEMS 3
Alternatively (1.1.3) can be rewritten as
x
i+1
= x
i
+x
i
(1.1.4)
where
x
i
= hf
i
+
h
2
2!
f
i
+
h
3
3!
f
i
+··· (1.1.5)
Starting from i = 0 with x
0
given and x
0
computed based on any desired
number of terms in (1.1.5), the value of x
1
is first calculated. Then, by letting
i = 1, 2, etc., in (1.1.4), the values of x
2
, x
3
, etc., are obtained successively. The-
oretically, if the number of terms retained in (1.1.5) increases indefinitely, the
numerical result from this marching scheme approaches the exact solution. In
reality, however, it is not permissible to do so, and the series has to be cut off
after a certain finite number of terms. For example, if two terms are retained on
the right-hand side of (1.1.5) in computing x
i
, the value of x
i+1
so obtained
is smaller than the exact value by an amount (h
3
3!)f
i
+(h
4
4!)f
i
+···. For
small h the first term is dominant. We may say that the error involved in this
numerical calculation is of the order of h
3
f
i
,orsimplyO(h
3
f
i
).Thisisthe
truncation error that results from taking a finite number of terms in an infinite
series.
It is termed Euler’s method when only one term is used on the right-hand side
of (1.1.5). The truncation error is O(h
2
f
i
), and the method should not be used if
accuracy is demanded in the result.
It is impractical to use Taylor’s series expansion method if f is a function that
has complicated derivatives. Furthermore, because of the dependence of the series
on the derivatives of f , a generalized computer program cannot be constructed for
this method. The nth-order Runge-Kutta method is a commonly used alternative.
Computations in this method require the evaluation of the function, f ,instead
of its derivatives, with properly chosen arguments; the accuracy is equivalent
to that with n terms retained in the series expansion (1.1.5). The second-order
Runge-Kutta formulas are
x
i+1
= x
i
+hf
x
i
+
1
2
1
x
i
, t
i
+
1
2
h
(1.1.6)
where
1
x
i
= hf (x
i
, t
i
) (1.1.7)
For better results the following fourth-order Runge-Kutta formulas are usually
employed:
x
i+1
= x
i
+
1
6
(
1
x
i
+2
2
x
i
+2
3
x
i
+
4
x
i
) (1.1.8)
in which the increments are computed in the following order:
1
x
i
= hf (x
i
, t
i
)
2
x
i
= hf
x
i
+
1
2
1
x
i
, t
i
+
1
2
h
3
x
i
= hf
x
i
+
1
2
2
x
i
, t
i
+
1
2
h
4
x
i
= hf (x
i
+
3
x
i
, t
i
+h)
(1.1.9)
4 FLOW TOPICS GOVERNED BY ORDINARY DIFFERENTIAL EQUATIONS
The derivation of these formulas can be found, for instance, in Kuo (1972,
p. 137). In the fourth-order method the dominant term in the truncation error
is (h
5
5!)f
iv
i
. Unless the slope of the solution is very steep, satisfactory results
are usually obtained with reasonably small h. Stability and step-size control in
Runge-Kutta methods are referred to in the book by Carnahan, Luther, and Wilkes
(1969, p. 363) and are not discussed here.
Runge-Kutta methods can be extended to solve higher-order or simultaneous
ordinary differential equations. Consider a second-order equation of the gen-
eral form
d
2
x
dt
2
= F
x,
dx
dt
, t
(1.1.10)
accompanied by the initial conditions that x = x
0
and dx
dt = p
0
when t = t
0
.
By calling the first-order derivative a new variable p, the equation (1.1.10) can
be written in the form of two simultaneous first-order equations:
dx
dt
= p
dp
dt
= F(x, p, t)
(1.1.11)
with initial values x(t
0
) = x
0
and p(t
0
) = p
0
. The fourth-order Runge-Kutta for-
mulas for marching from t
i
to t
i+1
are
x
i+1
= x
i
+
1
6
(
1
x
i
+2
2
x
i
+2
3
x
i
+
4
x
i
)
p
i+1
= p
i
+
1
6
(
1
p
i
+2
2
p
i
+2
3
p
i
+
4
p
i
)
(1.1.12)
in which the individual terms are computed in the following order:
1
x
i
= hp
i
1
p
i
= hF(x
i
, p
i
, t
i
)
2
x
i
= h
p
i
+
1
2
1
p
i
2
p
i
= hF
x
i
+
1
2
1
x
i
, p
i
+
1
2
1
p
i
, t
i
+
1
2
h
3
x
i
= h
p
i
+
1
2
2
p
i
3
p
i
= hF
x
i
+
1
2
2
x
i
, p
i
+
1
2
2
p
i
, t
i
+
1
2
h
4
x
i
= h(p
i
+
3
p
i
)
4
p
i
= hF(x
i
+
3
x
i
, p
i
+
3
p
i
, t
i
+h)
(1.1.13)
Most of the programs included in this chapter have utilized MATLAB initial-
value solvers, ODE23 and ODE45; however, we also provide explicit function
subprograms that contain fourth-order Runge-Kutta solvers. These subprograms
can be easily substituted for the MATLAB solvers by the user.
FREE FALLING OF A SPHERICAL BODY 5
1.2 FREE FALLING OF A SPHERICAL BODY
As the first application of the Runge-Kutta methods, the motion of a free-falling
body is to be studied. This problem is an example in which a solution cannot be
obtained without using a numerical method.
It is said that Galileo released simultaneously two objects of different masses
from the Leaning Tower of Pisa and found that they touched the ground at the
same instant. If Galileo did perform such an experiment, is the conclusion stated
in the story correct? Certainly it is true in a vacuum. In the atmosphere, however,
forces are exerted on a body by the surrounding air that are determined by the size
and motion of the body, and the conclusion seems doubtful. If it is not correct,
which body should touch the ground first, the larger one or the smaller one?
To answer these questions, we will not repeat the experiment in a laboratory.
Instead, we will first formulate the body motion, including the forces caused by
the surrounding fluid, and then perform the experiment numerically on a digital
computer. Because of its available drag data, a spherical body is preferred for
our analysis.
The z -axis is chosen in the direction of gravitational acceleration g; its origin
coincides with the center of the sphere at the initial instant t = 0, as shown in
Fig. 1.2.1. At time t > 0, the sphere of diameter d and mass m is at a distance z
from the origin and has a velocity v. It is surrounded by a fluid of density ρ
f
and
kinematic viscosity ν. In a vacuum the only external force acting on the body is
the gravitational pull, mg, in the positive z direction. While moving through a
fluid it is acted on by the following additional forces:
1. The buoyant force. According to Archimedes’ principle, the buoyant force
is equal to the weight of the fluid displaced by the body. It has the expression
−m
f
g,where
m
f
=
1
6
πd
3
ρ
f
(1.2.1)
t = 0
z = 0
dt
dz
v =
z
g
d
m
FIGURE 1.2.1 A free-falling spherical body.
6 FLOW TOPICS GOVERNED BY ORDINARY DIFFERENTIAL EQUATIONS
and the negative sign means that the force is in the direction of the negative
z-axis.
2. The force on an accelerating body. When a body immersed in a stationary
fluid is suddenly set in motion, a flow field is induced in the fluid. The kinetic
energy associated with the fluid motion is generated by “doing work,” or moving
the body against a drag force. This drag exists even if the fluid is frictionless,
and it has the value −
1
2
m
f
dv/dt, derived from an inviscid theory (Lamb, 1932,
p. 124).
3. The forces caused by viscosity. Around a body moving through a real fluid,
a region of rapid velocity change exists adjacent to the surface. The velocity
gradient causes a shear stress on the surface; the force resulting from integrating
the shear stresses throughout the body is termed the skin friction. In addition, in
a viscous fluid the pressure at the rear of the body becomes lower than that at the
front. The pressure difference gives rise to a pressure or form drag on the body.
The total viscous force, or the sum of skin friction and form drag, is sometimes
expressed in a dimensionless form called a drag coefficient, which is a function
of the body shape and the Reynolds number. Except for a few simple shapes
and at very low Reynolds numbers, this function is difficult to find analytically,
and it is usually determined through experiment. Figure 1.2.2 shows a typical
experimental curve for smooth spheres (Goldstein, 1938, p. 16) in which the
drag coefficient c
d
, defined as the total viscous force divided by
1
2
ρ
f
v
2
1
4
πd
2
,is
plotted against the Reynolds number Re, defined as v d
ν.
4. The wave drag. When the body speed is comparable to the speed of sound
in a fluid medium, shock waves may develop on or ahead of the body, causing
a wave drag.
FIGURE 1.2.2 Drag coefficient for smooth spheres.
FREE FALLING OF A SPHERICAL BODY 7
If we consider only low subsonic speeds, the wave drag can be safely omitted.
Including all the other forces discussed, Newton’s law of motion, when applied
to a spherical body, has the form
m
dv
dt
= mg −m
f
g −
1
2
m
f
dv
dt
−
1
2
ρ
f
v|v|
π
4
d
2
c
d
(v)
v|v| is used instead of v
2
, so that the direction of viscous drag is always opposite
to the direction of v. After rearranging, it can be written as
m +
1
2
m
f
dv
dt
= (m − m
f
)g −
π
8
ρ
f
v|v|d
2
c
d
(v) (1.2.2)
The left-hand side indicates that in an accelerating or decelerating motion
through a fluid, the body behaves as if its mass were increased. The term
1
2
m
f
is sometimes referred to as the added mass. Upon substitution from (1.2.1) and
from m =
1
6
πd
3
ρ, ρ being the density of the body, (1.2.2) becomes
dv
dt
=
1
A
[
B −C v|v|c
d
(v)
]
(1.2.3)
with
dz
dt
= v (1.2.4)
where A = 1 +
1
2
ρ, B = (1 − ρ)g,andC = 3ρ/4d; ρ stands for the density
ratio ρ
f
/ρ. These equations fit the general form (1.1.11) and can be solved by
applying Runge-Kutta methods.
Special care is needed for the empirical function c
d
(v) in numerical compu-
tation. To present it in a form suitable for the computer, the curve is replaced
by several broken lines, as shown in Fig. 1.2.2. To the left of the point a,where
Re ≤ 1, the Stokes formula
c
d
=
24
Re
(1.2.5)
is used, which coincides well with the experimental curve for Reynolds numbers
that are much less than unity and deviates only slightly from it in the neighbor-
hood of the point a (where Re = 1, c
d
= 24). On the log–log plot a straight line
is drawn between points a and b (where Re = 400, c
d
= 0.5) to approximate
the actual curve, which gives
c
d
= 24
Re
0.646
for 1 < Re ≤ 400 (1.2.6)
Between b and c (where Re = 3 ×10
5
), we assume that the drag coefficient
has a constant value of 0.5. The abrupt drop of drag coefficient around c is
an indication of transition from laminar to turbulent boundary layer before flow
separates from the body surface (see Kuethe and Chow, 1998, Section 17.11).
Another straight line is drawn between d (where Re = 3 ×10
5
, c
d
= 0.08) and
e (where Re = 2 × 10
6
, c
d
= 0.18). Thus,
c
d
= 0.000366 ×Re
0.4275
for 3 × 10
5
< Re ≤ 2 ×10
6
(1.2.7)
8 FLOW TOPICS GOVERNED BY ORDINARY DIFFERENTIAL EQUATIONS
Finally, in the high-Reynolds number region beyond e, a constant value of 0.18
is assumed for the drag coefficient.
Here we have chosen a simple way to approximate a complicated function.
Better results can be expected by dividing the curve into a larger number of
broken lines. Alternatively, the drag coefficient can be fed into the computer
as a tabulated function from which the value at any Reynolds number can be
interpolated or extrapolated.
Now let us return to the numerical experiment. Suppose steel spheres are
dropped in air under standard atmospheric conditions at sea level; we have
ρ = 8000 kg/m
3
, ρ
f
= 1.22 kg/m
3
, ν = 1.49 ×10
−5
m
2
/s, and g = 9.8m/s
2
.
Because of the low value of
ρ in this particular example, the effects of buoyancy
and added mass are negligible, as can be seen in (1.2.3). They become significant
if the experiment is repeated in water or in any fluid whose density is compara-
ble to that of the body. However, for generality these terms are still kept in our
computer program.
The result for a sphere 0.01 m in diameter will be shown. Its motion is to be
compared with that of the same body when dropped in a vacuum. In the latter
case
ρ = 0, (1.2.3) and (1.2.4) can be integrated directly to give the solution in
the vacuum:
v
v
= v
0
+gt
z
v
= z
0
+v
0
t +
1
2
gt
2
where z
0
and v
0
are, respectively, the initial position and velocity of the sphere at
t
0
= 0. We will choose the initial values z
0
= 0andv
0
= 0 in the computation.
Since the right-hand side of (1.2.3) does not explicitly contain z or t,letus
call it F (v), which has a form simpler than the general expression shown in
(1.1.11). We compute the position and velocity at time t
i+1
based on those at t
i
according to the following formulas simplified from (1.1.12) and (1.1.13).
The eight increments are first calculated.
1
z
i
= hv
i
1
v
i
= hF(v
i
)
2
z
i
= h
v
i
+
1
2
1
v
i
2
v
i
= hF
v
i
+
1
2
1
v
i
3
z
i
= h
v
i
+
1
2
2
v
i
3
v
i
= hF
v
i
+
1
2
2
v
i
4
z
i
= h(v
i
+
3
v
i
)
4
v
i
= hF(v
i
+
3
v
i
)
They enable us to compute the new values:
z
i+1
= z
i
+
1
6
(
1
z
i
+2
2
z
i
+2
3
z
i
+
4
z
i
)
v
i+1
= v
i
+
1
6
(
1
v
i
+2
2
v
i
+2
3
v
i
+
4
v
i
)
FREE FALLING OF A SPHERICAL BODY 9
According to the above formulas, the variables t
i
, z
i
,andv
i
are one-
dimensional arrays whose elements consist of their individual values at
t
0
, t
1
, t
2
, .... However, if instead of printing the complete set of data at the last
time step we print out the time, position, and velocity at every time step t
i
,the
same variable names may be used for the quantities evaluated at t
i+1
.Inthis
way the subscripts i and i +1 can be dropped from the variables to simplify
the appearance of the program.
In Program 1.1 the instantaneous Reynolds number is shown in addition to
the position and velocity for our reference. A time limit
TMAX is introduced in
the program. The numerical integration is stopped when t exceeds this value.
In Program 1.1 a time step of 0.1 s is used. To check the accuracy, the same
program has been repeated with step sizes of 0.02 and 0.2 s. Since no appreciable
difference can be found between the three sets of data, it can be concluded that
the result is reliable. Keeping h = 0.1 s, more data are obtained by varying the
diameter while doubling the maximum time of integration. They are plotted in
Figs. 1.2.3 and 1.2.4.
We have thus obtained some results from the numerical experiment of a free-
falling steel sphere. The displacement curves in Fig. 1.2.3 indicate that to travel
through a given distance, a larger body takes less time than a smaller one. For
two spheres of comparable sizes, the difference in arrival times at a distance of
56 m (the height of the Tower of Pisa) is only a small fraction of a second. Such
a time difference is difficult to detect without the help of some instruments.
0 2 4 6 8 10 12
0
50
100
150
200
250
300
350
400
450
500
t (s)
z (m)
d = 0.001 m
0.005
0.01
0.02
0.07
In vacuum
or d → ∞
FIGURE 1.2.3 Displacement of steel spheres falling in air.
10 FLOW TOPICS GOVERNED BY ORDINARY DIFFERENTIAL EQUATIONS
0 2 4 6 8 10 12
0
10
20
30
40
50
60
70
80
90
100
t (s)
v (m/s)
d = 0.001 m
0.005
0.01
0.02
0.07
In vacuum
or d → ∞
FIGURE 1.2.4 Velocity of steel spheres falling in air.
A body falling through a fluid always reaches a constant velocity, called the
terminal velocity, which increases with the diameter of a spherical body, as
shown in Fig. 1.2.4. For a small sphere the Reynolds number is relatively low
and the boundary layer always remains laminar before the point of separation.
At a sufficiently high Reynolds number, when the transition from a laminar to
a turbulent boundary layer occurs on a larger sphere, the abrupt decrease in
drag shown in Fig. 1.2.2 will cause a sudden acceleration of the body. This phe-
nomenon can be observed in Fig. 1.2.4 on the curve for a sphere with d = 0.07 m
around t = 7.5 s. For an extremely large sphere the effect of the surrounding fluid
becomes negligible in comparison with body inertia, so that the sphere behaves
as if it were moving in a vacuum. In this case the velocity increases indefinitely
with time, and a terminal velocity can never be reached.
When traveling at the terminal velocity, the gravitational force is balanced by
the sum of buoyancy and viscous drag, so that the acceleration is zero. From
(1.2.2) an expression is derived for the terminal velocity:
v
t
=
4(1 −ρ)gd
3ρc
d
(v
t
)
(1.2.8)
Based on the numerical result for the steel sphere of 0.01 m diameter, the terminal
Reynolds number is estimated to be not much higher than 2.73 × 10
4
. According
FREE FALLING OF A SPHERICAL BODY 11
to the approximation shown in Fig. 1.2.2, the drag coefficient at the terminal
velocity is c
d
(v
t
) = 0.5. Substitution of this together with other given values
into (1.2.8) gives v
t
= 41.4m/s.
In general, c
d
(v
t
) is not known a priori, and the terminal velocity cannot be
obtained unless the time history of the motion is computed first. However, from
the data computed for steel spheres falling through air, it is concluded that the
terminal Reynolds number exceeds 2 ×10
6
if the diameter of a sphere is greater
than 0.12 m. c
d
(v
t
) has a constant value of 0.18 for such a sphere according to
our approximation, and v
t
can readily be calculated. On the other hand, for tiny
particles whose terminal Reynolds numbers are in the Stokes flow regime, the
Stokes formula (1.2.5) is used in (1.2.8) to obtain
v
t
=
(1 −
ρ)gd
2
18ρν
for Re < 1 (1.2.9)
Let us now estimate the minimum diameter of a steel sphere whose terminal
velocity in air can exceed the speed of sound, 340 m/s at sea level. The terminal
Reynolds number of such a body is well beyond 2 ×10
6
, so that the value 0.18
is used for c
d
(v
t
). By assigning the sound speed to v
t
, (1.2.8) gives d = 0.243 m.
Remember that before reaching sonic speed, a transonic region appears on the
body, causing a wave drag that should also be included in the drag coefficient.
Furthermore, because of the limited drag data available for high Reynolds num-
bers, the constant c
d
assumption is doubtful. The minimum diameter estimated
here is expected to be too low.
In the atmosphere or in an ocean, if the displacement of a body is so large
that the variations in fluid density and kinematic viscosity are significant, the
program must be modified by specifying the dependence of these quantities on
the height z. In this case the function F(v) in Program 1.1 should be changed to
F(v, z).
Problem 1.1 Find the motion of a Ping-Pong ball 0.036 m in diameter released
at the bottom of a water tank. The density of water is 1000 kg/m
3
and the
kinematic viscosity is 1 × 10
−6
m
2
/s. Assume that the shell of the ball is so thin
that the density of the body is the same as that of the filling air, or 1.22 kg/m
3
.
Hint
Because of the large buoyant force, the stationary ping-pong ball initially expe-
riences a large acceleration. You may verify that the numerical solution diverges
if h = 0.1 s is still used. The value h = 0.01 s is recommended in this problem.
Since the ball reaches an upward steady motion within a fraction of a second, a
maximum time of 0.5 s is sufficient for the computation.
Problem 1.2 For a given body falling in a specified fluid, the terminal velocity
is fixed independent of the initial velocity of the body. If the initial velocity is
slower than or in a direction opposite to the terminal velocity, the body will accel-
erate toward the terminal value. On the other hand, if the initial velocity is faster